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Abstract 


The  ability  to  produce  high-resolution  images  of  the  Earth’s  surface  from  space 
has  flourished  in  recent  years  with  the  continuous  development  and  improvement 
of  satellite-based  imaging  sensors.  Earth-imaging  satellites  often  rely  on  complex 
onboard  navigation  systems,  with  dependence  on  Global  Positioning  System  (GPS) 
tracking  and/or  continuous  post-capture  georegistration,  to  accurately  geolocate  ground 
targets  of  interest  to  either  commercial  and  military  customers.  Consequently,  these 
satellite  systems  are  often  massive,  expensive,  and  susceptible  to  poor  or  unavail¬ 
able  target  tracking  capabilities  in  GPS-denied  environments.  Previous  research  has 
demonstrated  that  a  tightly-coupled  image-aided  inertial  navigation  system  (INS), 
using  existing  onboard  imaging  sensors,  can  provide  significant  target  tracking  im¬ 
provement  over  that  of  conventional  navigation  and  tracking  systems.  Satellite-based 
image-aided  navigation  is  explored  as  a  means  of  autonomously  tracking  stationary 
ground  targets  by  implementing  feature  detection  and  recognition  algorithms  to  accu¬ 
rately  predict  a  ground  target’s  pixel  location  within  subsequent  satellite  images.  The 
development  of  a  robust  satellite-based  image-aided  INS  model  offers  a  convenient, 
low-cost,  low- weight  and  highly  accurate  solution  to  the  geolocation  precision  prob¬ 
lem,  without  the  need  of  human  interaction  or  GPS  dependency,  while  simultaneously 
providing  redundant  and  sustainable  satellite  navigation  capabilities. 
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I.  Introduction 

This  research  outlines  the  methodology  of  integrating  the  imaging  sensors  of  a 
satellite-based  Earth  observation  system  with  the  inertial  sensors  of  the  satel¬ 
lite’s  navigation  system  in  order  to  accurately  locate  a  ground  target  of  interest  and 
autonomously  track  that  target  over  an  extended  period  of  time.  The  image-aided 
satellite  navigation  system  design  uses  conventional  feature  detection  and  recognition 
methods  in  order  to  provide  robust,  rapid  target  tracking  capability  without  the  need 
for  additional  vehicle-based  or  ground-based  hardware  or  dependency  on  external 
navigation  reference  sources. 

1.1  Background 

As  early  as  the  mid-1950s,  the  desire  to  explore  the  space  beyond  our  skies  led 
to  the  development  of  man-made  satellites.  Ever  since  the  first  of  these  satellites,  the 
Soviet-built  Sputnik  1,  was  launched  in  1957,  scientists,  astronomers,  and  engineers 
were  eager  to  harness  the  vast  capabilities  of  the  space  environment.  These  efforts 
paved  the  way  for  space-based  breakthroughs  in  communications,  weather,  imaging, 
and  manned  space  flight. 

In  particular,  the  numerous  applications  involving  terrestrial  imagery  from  space 
have  led  to  high  demands  in  the  areas  of  agriculture,  geology,  education,  intelligence, 
and  warfare.  At  any  one  point  in  time,  day  or  night,  ground  images  are  continuously 
captured  from  space  by  observation  satellite  systems  orbiting  the  Earth.  Often,  a 
ground  target  of  interest’s  geographic  location  (referred  to  as  the  target’s  “geoloca- 
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tion”)  is  desired;  for  example,  for  commercial  or  military  purposes.  However,  the 
ability  to  quickly,  accurately  and  efficiently  detect  and  track  a  ground  target  from 
space,  without  the  need  of  sophisticated,  costly,  or  deniable  systems,  has  been  an 
ongoing  challenge  of  the  space  community. 

1.2  Problem  Definition 

Current  Earth-imaging  satellite  systems  often  rely  on  either  external  navigation 
reference  sources,  such  as  the  Global  Positioning  System  (GPS),  or  continuous  post- 
capture  georegistration,  to  provide  precise  ground  target  geolocation.  Often,  these 
current  imaging  system  designs  are  massive,  expensive,  sluggish,  and  jammable.  An 
additional  risk  to  highly  complex  systems  is  the  increased  susceptibility  to  unforeseen 
hardware-related  anomalies.  Since  high  image  resolution  requirements  must  be  met  in 
order  to  track  many  man-made  targets,  conventional  satellite  databases  often  demand 
high  memory  capacity.  The  resulting  image  processing  time  is  often  impractical. 
Additionally,  if  routine  hardware  calibration  is  required  onboard  the  satellite  in  order 
to  maintain  navigation  accuracy,  the  resulting  image  capturing  performance  could 
degrade  if  such  a  calibration  where  to  inhibit  the  imaging  sensor’s  ability  to  locate 
the  ground  target. 

The  motivation  of  this  research  is  to  develop  an  inexpensive,  light-weight,  highly 
accurate  image-aided  inertial  satellite  navigation  system  without  the  need  for  human 
interaction  or  dependence  on  external  navigation  reference  systems,  such  as  GPS. 
Image-aided  geolocation  algorithms  are  of  interest  for  both  military  and  civilian  space 
applications  in  which  weight,  cost,  and  time  savings  are  of  interest  and/or  the  denial  of 
GPS  is  of  concern.  Any  civilian  or  military  satellite  system  producing  digital  imagery 
would  benefit  from  a  robust,  low-cost,  autonomous  satellite  system  providing  precise 
geolocation  capabilities. 
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1 . 3  Scope 

Satellite-based  image-aided  navigation  is  explored  as  a  means  of  autonomously 
tracking  a  stationary  ground  target  within  subsequent  satellite  images  and  using  fea¬ 
ture  recognition  algorithms  to  predict  the  target’s  pixel  location  within  each  image. 
In  order  to  make  this  problem  tractable,  certain  assumptions  are  made.  Although  the 
effects  of  atmospheric  turbulence  are  covered  in  detail,  weather  and  lighting  conditions 
are  assumed  to  be  favorable  at  the  location  of  interest.  For  example,  environmental 
limitations  such  as  precipitation,  fog,  cloud  cover,  and  poor  sunlight  conditions  are 
not  covered.  Image  quality  issues,  including  ego-motion  disparity  and  motion  blur 
among  subsequent  images  are  also  not  specifically  identified  in  this  research.  Also, 
for  simplicity,  a  spherical  Earth  model  of  constant  elevation  is  assumed. 

1.4  Research  Contributions 

Many  contributions  are  made  on  behalf  of  this  research.  These  contributions 
are  primarily  defined  in  Chapter  II  and  implemented  in  Chapter  III  of  this  thesis,  and 
are  listed  here  in  their  respective  order  of  appearance. 

The  extended  Kalman  filter  is  utilized  to  efficiently  estimate  the  state  of  the  non¬ 
linear  dynamic  satellite  navigation  system  from  a  series  of  noisy  measurements  [12,13]. 
A  background  in  orbital  mechanics,  specifically,  the  systematic  transformation  from 
geometrical  orbital  parameters  to  satellite  position  and  velocity  vectors,  provide  nom¬ 
inal  satellite  trajectory  modeling  [15,22],  A  detailed  understanding  of  the  turbu¬ 
lent  and  refractive  effects  on  an  image  propagating  through  the  atmosphere  provides 
insight  to  realistic  image  system  modeling  [1,5,21],  The  primary  contribution  to 
image-aided  navigation  theory  demonstrates  that  tightly-coupled  image-aided  inertial 
navigation  can  provide  significant  target  tracking  improvement  over  that  of  conven¬ 
tional  navigation  systems  [28,30],  Image  matching  theory  utilized  the  sum-of-squared- 
difference  algorithm,  in  which  a  target  with  unique  features  can  be  tracked  among  a 
series  of  pixclated  images  by  measuring  the  correspondence  between  images  [3,16], 
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Finally,  the  development  of  the  stochastic  projection  method  provided  a  means  of  ac¬ 
curately  and  optimally  predicting  a  target’s  pixel  location  within  in  a  series  of  images 
by  limiting  the  number  false  matches  and  constraining  the  feature  correspondence 
search  [28,30]. 

1.5  Methodology 

The  thesis  is  organized  as  follows: 

•  Chapter  II:  Chapter  II  provides  the  detailed  mathematical  background  of  the 
image/inertial  sensor  integration  problem.  The  first  sections  of  this  chapter 
cover  mathematical  notation,  reference-frame  declarations,  coordinate  transfor¬ 
mations,  inertial  sensor  design,  and  Kalman  filtering.  The  following  sections 
review  orbital  dynamics,  specifically,  Newtonian  and  Keplerian  theory,  satellite 
trajectory  transformation,  and  orbital  classifications.  Further  sections  discuss 
optical  modeling,  spatial  coherence,  and  imaging  through  turbulence,  specifi¬ 
cally,  the  effects  of  image  jitter  and  horizontal  light  refraction  displacement. 
The  following  sections  cover  image-aided  navigation  theory,  image  matching 
techniques,  and  georeferencing  theory.  The  final  section  of  this  chapter  covers 
the  stochastic  projection  method. 

•  Chapter  III:  Chapter  III  covers  the  methodology  of  this  research,  implementing 
the  theories  covered  in  Chapter  II  to  build  a  satellite-based  image-aided  nav¬ 
igation  system.  The  first  sections  develop  the  orbital  model,  satellite  vehicle 
model,  and  imaging  system  model.  The  next  sections  assign  noise  modeling 
with  respect  to  atmospheric  turbulence,  image  sensors,  satellite  trajectory,  and 
optical  measurements.  The  following  sections  develop  the  truth  model  with  re¬ 
spect  to  the  ground  target  and  implement  image  matching  and  georeferencing 
techniques  to  provide  an  accurate  target  location  error  prediction.  The  final 
section  of  this  chapter  is  the  implementation  of  the  extended  Kalman  filter. 
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•  Chapter  IV :  Chapter  IV  supports  the  Monte  Carlo  results  and  observations  of 
the  image-aided  satellite  navigation  system.  Two  distinct  profiles,  one  satellite 
in  a  low-Earth  orbit  with  high  image  resolution,  and  the  other  in  a  high-Earth 
orbit  with  low  image  resolution,  are  produced.  The  respective  vehicle  position, 
vehicle  velocity,  vehicle  attitude,  and  target  location  errors  are  analyzed  in  detail 
in  the  following  sections,  both  with  and  without  the  introduction  of  image-aided 
target  predictions. 

•  Chapter  V :  Chapter  V  provides  conclusions  and  closing  remarks  regarding  the 
image/inertial  integrated  navigation  system,  as  well  as  potential  areas  for  future 
exploration  in  the  subject. 
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II.  Background 


This  chapter  reviews  the  mathematical  and  conceptual  background  required  to 
fully  develop  an  image-aided  navigation  system  of  an  orbiting  satellite.  First, 
a  definition  of  the  mathematical  notation  used  throughout  the  document  will  be 
presented,  followed  by  reference  frame  definitions.  A  basic  understanding  of  inertial 
navigation  will  follow.  Next,  a  review  of  both  linear  and  nonlinear  Kalman  filtering 
methods  will  be  discussed,  specifically,  the  Extended  Kalman  Filter  (EKF).  Orbital 
mechanics  required  to  define  a  satellite’s  orbital  path,  as  well  as  the  concepts  behind 
satellite  imaging  will  be  described.  Finally,  an  in-depth  discussion  of  imaging  through 
atmospheric  turbulence  will  be  presented  and  analyzed. 

2.1  Mathematical  Notation 

The  mathematical  notation  to  be  used  throughout  this  paper  is  listed  in  Ta¬ 
ble  2.1. 

2.2  Inertial  Navigation 

In  this  section,  basic  concepts  of  an  inertial  navigation  system  (INS)  are  dis¬ 
cussed,  including  navigation  reference  frames,  coordinate  frame  transformation,  and 
functionality  and  errors  associated  with  strapdown  INS  sensors. 

2.2.1  Basic  Concepts  .  The  concept  of  navigation  as  a  means  of  determining 
direction  from  one  place  to  another  has  been  used  for  centuries.  Navigation  can  be 
as  simple  as  following  directions  on  a  map  by  determining  position  based  on  one’s 
surroundings.  Navigation  systems  are  often  developed  for  vehicles  in  order  to  plot, 
ascertain,  and  direct  the  vehicle  through  land,  air,  sea  or  space.  One  navigation 
technique  uses  fixed  stars  to  define  a  reference  frame  fixed  in  space.  This  reference 
is  commonly  referred  to  as  the  “inertial”  reference  frame.  Given  knowledge  of  the 
motion  of  the  Earth  and  the  time  of  the  observation,  the  navigator  is  able  to  use  the 
celestial  measurements  to  define  his  or  her  position  on  the  surface  of  the  Earth  [24], 
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Table  2.1:  Defined  mathematical  notation 


Type 

Description 

Example 

Scalars 

Scalar  variables  are  designated  with  italic  type 

x  or  X 

Vectors 

Vectors  are  denoted  by  lower  case  bold  type 

X 

Matrices 

Matrices  are  denoted  by  uppercase  bold  type 

X 

Transpose 

The  transpose  of  a  vector  or  matrix  is  desig¬ 
nated  with  a  superscript  capital  letter  T 

x2  or  X7 

Estimated  Variables 

Estimates  of  random  variables  are  identified 
with  the  hat  character 

X 

Calculated  Variables 

Variables  containing  error  are  denoted  with 
the  tilde  character 

X 

Direction  Cosine  Ma¬ 
trix  (DCM) 

DCMs  are  designated  by  a  bold  capital  letter 
C  with  a  subscript  designating  the  originat¬ 
ing  coordinate  frame  and  a  superscript  as  the 
resulting  coordinate  frame 

c* 

Frame  of  Reference 

Vectors  expressed  in  a  specific  reference  frame 
are  annotated  with  a  superscript  letter 

j.a 

2.2.2  Reference  Frames.  In  order  to  express  inertial  navigation  information 
in  standardized  coordinates,  fundamental  reference  frames  must  be  defined  relative 
to  an  origin  and  orthogonal  axes.  Descriptions  of  the  reference  frames  discussed  are 
as  follows,  based  on  those  described  in  Refs.  [24]  and  [25]: 

•  The  True  Inertial  Frame  (/-frame)  is  a  theoretical  reference  frame  where  New¬ 
ton’s  laws  of  motion  apply;  therefore,  it  has  no  predefined  origin  or  orientation. 

•  The  Earth- centered  Inertial  Frame  (i-frame),  depicted  in  Figure  2.1,  is  an  or¬ 
thogonal  reference  frame  with  an  origin  at  the  Earth’s  center  of  mass  and  a 
non-rotating  x,  y  and  z  axes  with  respect  to  fixed  stars.  For  terrestrial  naviga¬ 
tion  purposes,  the  Earth-centered  inertial  frame  can  be  considered  an  inertial 
reference  frame. 

•  The  Earth- centered  Earth-fixed  Frame  (e-frame),  depicted  in  Figure  2.2,  is  an 
orthogonal  reference  frame  whose  origin  is  also  located  at  Earth’s  center  of 
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mass.  Its  x,  y  and  z  axes  are  fixed  with  respect  to  the  Earth,  with  the  x  axis 
on  the  equatorial  plane  pointing  toward  the  Greenwich  meridian,  the  y  axis  on 
the  equatorial  plane  pointing  toward  90  degrees  east  longitude,  and  the  z  axis 
aligned  with  the  north  pole. 

•  The  Vehicle-fixed  Navigation  Frame  (n'- frame) ,  depicted  in  Figure  2.1,  is  an 
orthogonal  reference  frame  with  an  origin  located  at  a  predefined  point  on  a 
vehicle  (usually  the  vehicle’s  center  of  gravity).  Its  x,  y  and  z  axes  point  in  the 
north,  east  and  down  (NED)  directions  relative  to  the  Earth,  respectively.  For 
standardization,  down  is  defined  as  the  direction  of  the  local  vertical  component 
of  the  Earth’s  gravity  vector. 

•  The  Earth-fixed  Navigation  Frame  (n-frame),  depicted  in  Figure  2.2,  is  an  or¬ 
thogonal  reference  frame  with  an  origin  located  at  a  predefined  location  on  the 
Earth  (i.e. ,  Earth’s  surface).  Similar  to  the  n ’-frame,  its  x,  y  and  z  axes  point 
in  the  NED  directions  relative  to  the  Earth,  respectively  (where  down  is  defined 
as  the  direction  of  the  local  vertical  component  of  Earth’s  gravity  vector). 

•  The  Body  Frame  (b- frame),  depicted  in  Figure  2.3,  is  an  orthogonal  reference 
frame  aligned  with  the  roll,  pitch  and  yaw  axes  that  point  out  of  a  vehicle’s 
nose,  right  wing  and  bottom,  respectively.  Vehicle  strapdown  inertial  navigation 
systems  are  referenced  in  the  6-frame.  For  an  orbiting  satellite,  it  is  assumed 
that  this  frame  is  equivalent  to  the  n-frame. 

•  The  Camera  Frame  (c- frame),  depicted  in  Figure  2.4,  is  an  orthogonal  reference 
frame  rigidly  attached  to  a  camera,  with  origin  at  the  camera’s  optical  center.  Its 
x,  y  and  z  axes  are  oriented  up,  to  the  right,  and  out  of  the  camera,  respectively. 
It  will  be  assumed  in  this  research  that  the  camera  is  rigidly  mounted  onto  the 
satellite;  therefore,  the  c-frame  will  be  equivalent  to  the  n-frame. 

2.2.3  Coordinate  Transformations.  Coordinate  transformations  describe  the 
relationship  between  two  reference  frames  and  are  classified  as  either  three-  or  four- 
parameter  transformations  [19,24,25].  In  this  research,  this  coordinate  transformation 


Local 

Meridian 


ys 


Figure  2.1:  Inertial,  Earth  and  Navigation  Reference  Frames  [25]. 


Figure  2.2:  Earth-centered,  Earth-fixed  Navigation  Reference  Frame  [25]. 
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Figure  2.3:  Body  Reference  Frame  [24],  Figure  2.4:  Camera  Reference  Frame  [25]. 


will  be  the  direction  cosine  matrix  (DCM).  Since  three-parameter  coordinate  trans¬ 
formations  contain  a  singularity  at  a  pitch  angle  of  90°,  the  use  of  the  four-parameter 
DCM  coordinate  transformation  will  be  used. 

The  DCM  is  a  three-by-three  matrix  representing  the  unit  vector  of  the  origi¬ 
nating  frame  projected  along  the  axis  of  the  resulting  frame.  The  DCM  is  written  in 
component  form  as 


C 11  C12  C13 


(2.1) 


o  C21  C22  C23 


C31  C32  C33 


The  elements  ctJ  represent  the  cosine  of  the  angle  between  the  z-axis  of  originating 
reference  frame  (o)  and  the  j-axis  of  the  resulting  frame  (r).  The  DCM  is  propagated 
in  time  through  the  equation 


(2.2) 


where  fl°ro  is  the  skew  symmetric  form  of  the  angular  rate  vector  u°0  =  [ux  u)y  ujz]T, 
defined  as 
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(2.3) 


f! 


o 

ro 


0  —  LUZ  Uy 

L Oz  0  ~0JX 


(dy  OJ x  0 


This  represents  the  angular  turn  rate  of  the  originating  frame  with  respect  to  the 
resulting  frame  expressed  in  the  axes  of  the  originating  frame. 


2.2.4  Navigation  Sensors  .  To  more  accurately  describe  inertial  navigation, 
it  is  the  process  of  establishing  the  position,  velocity,  and  attitude  of  a  vehicle  using 
information  derived  from  its  inertial  sensors.  Sensors  within  INS  units  are  often  rigidly 
“fixed”  to  the  moving  body,  commonly  referred  to  as  a  strapdown  INS  system.  These 
sensor  devices,  namely  accelerometers  and  gyroscopes,  are  used  to  measure  linear  and 
angular  motion,  respectively,  in  the  inertial  frame.  Figures  2.5  and  2.6  illustrate  the 
function  of  accelerometers  and  gyroscopes  within  a  strapdown  INS  system. 

Accelerometers  are  mechanical  or  electrical  sensor  systems  that  use  seismic 
masses  and  springs  to  measure  translational  motion  of  the  platform  in  which  they 
are  located.  In  other  words,  accelerometers  use  simple  proof  mass  physics  to  measure 
the  total  external  specific  force  acting  upon  itself  [24],  Specific  force  is  defined  here 
as  the  sum  of  acceleration  acting  on  the  body  with  respect  to  the  inertial  reference 
frame  plus  gravity.  Strapdown  IMU  systems  consist  of  a  triad  of  accelerometers, 
usually  aligned  with  the  vehicle’s  body  reference  frame. 

Gyroscopes  are  the  mechanical,  electrical,  or  optical  sensor  systems  used  to 
measure  rotational  motion  in  an  inertial  system  [24],  Mechanical  gyroscopes  rely  on 
the  inertial  properties  of  a  proof  (often  spinning)  mass  for  their  operation,  producing 
measurements  of  turn  angle  or  turn  rate  with  respect  to  inertial  space.  Similarly, 
optical  gyroscopes  provide  a  measure  of  angular  rate.  However,  instead  of  a  using  a 
spinning  mass  to  detect  rotation,  these  gyroscopes  use  interference  of  laser  light,  called 
the  Sagnac  effect,  to  detect  changes  in  orientation  and  spin  of  the  gyroscope.  Optical 


11 


Figure  2.5:  A  two-dimensional  strapdown  INS  unit  [24],  The  body  attitude,  6 ,  is 
computed  by  integrating  the  measured  angular  rate,  ojyb,  and  is  used  to 
resolve  the  specific  forces,  fxb  and  fzb,  into  the  reference  frame. 

gyroscopes,  such  as  ring  laser  gyroscopes,  have  an  advantage  over  their  mechanical 
counterparts,  as  there  are  no  moving  parts,  no  inherent  drift  due  to  friction,  and 
generally  compact  in  size  and  lightweight. 

As  with  all  real-world  IMUs,  the  measurements  from  accelerometers  and  gyro¬ 
scopes  are  corrupted  by  errors.  Most  of  these  errors  are  correctable  through  factory 
calibration  techniques;  however,  it  is  not  possible  to  remove  all  errors.  A  brief  ex¬ 
planation  of  these  sensor  errors,  general  to  both  accelerometers  and  gyroscopes,  are 
listed  below  [24]  [25]: 

•  Bias'.  Constant  or  slowly- varying  additive  error. 

•  Scale  Factor :  Constant  or  slowly-varying  multiplicative  error. 

•  Sensor  Misalignment :  The  result  of  mechanical  fabrication  and  installation  er¬ 
rors  made  at  the  factory.  These  errors  result  in  a  difference  between  the  sensor’s 
sensitive  axis  and  the  platform  reference. 

•  Vibration :  Measurement  bias  as  a  function  of  specific  vibration(s). 
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Key: 


Figure  2.6:  A  strapdown  INS  unit  in  a  rotating  reference  frame  [24],  An  estimate 
of  the  turn  rate  of  the  reference  frame  is  derived  using  the  estimated 
component  of  horizontal  velocity. 


•  Measurement  Noise:  An  additive  error  component  with  high-bandwidth  power 
spectral  density,  such  as  electrical  noise,  thermal  noise,  etc. 


2. 3  Kalman  Filtering 

In  this  section,  basic  concepts  of  Kalman  filtering  are  discussed.  Two  types 
of  Kalman  filters  will  be  presented:  conventional  (for  linear  navigation  states)  and 
extended  (for  nonlinear  navigation  states).  All  equations  in  Section  2.3  are  derived 
in  Refs.  [12]  and  [13]. 

2.3.1  The  Conventional  Kalman  Filter.  The  conventional  Kalman  filter  is 
an  efficient  recursive  filter  that  estimates  the  state  of  a  dynamic  system  from  a  se¬ 
ries  of  measurements  containing  random  noise  [12],  To  make  the  estimation  problem 
tractable,  it  is  assumed  that  the  prior  knowledge  of  the  navigation  state  can  be  ade¬ 
quately  described  as  a  multivariate  Gaussian  distribution.  In  addition,  the  stochastic 
process  noise  and  additive  measurement  noise  are  assumed  to  be  zero-mean,  Gaussian 
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and  white,  and  that  the  nonlinear  state  dynamics  and  measurement  equations  can  be 
modeled  using  perturbation  techniques. 

The  Kalman  filter  state  equation,  also  known  as  the  linear  stochastic  differential 
equation,  can  be  written  as 

±(t)  =  F(t)x(t)  +  B(t)u(t)  +  G(t)w(t),  (2.4) 

where  F(t)  is  the  state  model  applied  to  the  state  vector  x(4),  B(t)  is  the  control- 
input  model  applied  to  the  control  vector  u(t),  and  G  (t)  is  the  noise  model  applied 
to  the  zero  mean,  white  Gaussian  process  noise  vector  w (t)  with  covariance  kernel 

E[ w(t)w7  (t  +  r)]  =  Q (t)h(r).  (2.5) 

In  Equation  (2.5),  Q (t)  is  the  process  noise  intensity  and  h(r)  represents  the  Dirac 
delta  function.  The  measurement  model  equation  for  the  Kalman  filter  is  defined  at 
time  tk  as 


z  (4)  =  H(4)x(4)  +  v(4),  (2.6) 

where  Ff(4)  is  the  observation  model  applied  to  the  true  state  vector  x(4)  and  v(4) 
is  the  zero-mean,  white  Gaussian  measurement  noise  process  of  intensity  R (i*)  in 
which 


F[v(^)vt(4)] 


J  R (ti)  U  =  tj 

1^0  ti  ^  tj. 


(2.7) 


The  initial  conditions  of  this  Kalman  filter  are  characterized  by  the  equations 


x(4)  =  x0 


(2.8) 
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and 


p  {to)  =  Po  (2.9) 

In  Equations  (2.8)  and  (2.9),  x0  is  defined  as  the  initial  measurement  estimate  vector 
at  initial  time  t0  and  PG  is  the  initial  covariance  matrix  of  the  form 

P  0  =  E[8xo6xZ],  (2.10) 

where  6xa  is  the  initial  state  error.  The  measurement  estimate  and  covariance  are 
propagated  forward  in  time  by  k  iterations  using  the  equations 


x(tfc+1)  =  <h(4+i,4)x(f+)  +  Bd(tk)n(tk)  (2.11) 

P(^fc+i)  =  ^(^fc+15  ^fc)P(^fc  )<^>t(4'+i,  tf.)  +  Gd(tk)Qd(tk)Gj(tk).  (2.12) 

In  Equations  (2.11)  and  (2.12),  x(tk+1)  and  P(tk+1)  are  the  propagated  state  vector 
and  covariance  matrix  prior  to  a  measurement  update,  respectively,  $(4+i,4)  = 
eF0fc)Ai  jg  the  state  transition  matrix,  B d(tk)  is  the  discrete  control-input  model, 
Gd(tk)  is  the  discrete  noise  model  and  Q  (tk)  is  the  discrete  process  noise  intensity 
(calculated  using  the  Van  Loan  approach  or  first-order  approximation  methods  [12]). 
The  measurement  updates  are  computed  as 


K(4)  =  P(t*)HT(tfc)  [H(4)P(4T)Hr(4)  +  R(4)]_1  (2.13) 

*(ffc)  =  x(tfc)  +  K(4)  [zfc  -  H(4)x(t")]  (2.14) 

P(tf )  =  P(^)  -  K(4)H(4)P(4),  (2.15) 
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where  K (£*.)  is  the  Kalman  gain,  zk  is  the  measurement,  and  x(t£)  and  P (tk)  are  the 
state  vector  and  covariance  matrix  following  the  measurement  update. 

2.3.2  The  Extended  Kalman  Filter  .  The  extended  Kalman  filter  (EKF)  is 
a  minimum  mean-square-error  (MMSE)  estimator  based  upon  the  first-order  approx¬ 
imation  Taylor  series  expansion  of  the  non-linear  system  dynamics  and  measurement 
models  [13].  Unlike  the  conventional  Kalman  filter,  the  EKF  can  handle  nonlinear 
system  models  and  is  therefore  more  applicable  in  real-world  scenarios.  The  EKF 
nonlinear  stochastic  differential  equation  can  be  written  as 

x(t)  =  f  [  x(t),  u(t),t  ]  +  G(t)w(t),  (2.16) 

where  f  [  x(t),u(t),t  ]  is  a  known  model  vector  of  nonlinear  functions  with  respect 
to  the  state  vector  x(t)  and  the  control-input  vector  u(t).  G (t)  is  the  noise  model 
applied  to  the  additive  noise  vector  w (t). 

Rewriting  the  nonlinear  measurement  equations  in  matrix  form  with  respect 
to  the  desired  navigation  states,  the  discrete-time  measurements  for  the  EKF  are 
modeled  as  a  known  nonlinear  function  of  the  state  plus  linearly  additive  measurement 
noise  as 


z (ti)  =  h  [  x(ti),ti  ]  +  v(U).  (2.17) 

Given  the  nonlinear  stochastic  state  and  measurement  models,  the  extended 
Kalman  filter  can  be  built.  As  a  starting  point,  the  initial  state  and  covariance 
conditions  at  time  to  are  defined  with  respect  to  the  initial  time  to-  They  are  defined 
as,  respectively, 


±(t„/t„)  =  x(t+) 


(2.18) 


and 
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P(«„/«o)  =  P((,+). 


(2.19) 


The  basic  concept  of  the  EKF  is  to  relinearize  about  each  state  estimate  x  once 
it  has  been  computed.  As  soon  as  a  new  state  estimate  is  made,  a  new  and  statistically 
“better”  reference  state  trajectory  is  incorporated  into  the  estimation  process.  With 
the  EKF,  it  is  assumed  that  deviations  from  the  reference  (or  nominal)  trajectory  are 
small  enough  to  allow  linear  perturbation  techniques  to  be  employed  with  adequate 
results. 

To  achieve  the  final  form  of  the  extended  Kalman  filter,  the  most  recent  nominal 
xn(f/fj),  defined  at  time  t  with  respect  to  incremented  time  £*,  is  combined  with  the 
state  perturbation  estimate,  5it(t/ti),  to  generate  an  estimate  of  the  full  state.  This 
assumes  the  following  model  is  accurate  for  x(f): 

x(t)  =  xn(f/fj)  +  Sx(t/ti).  (2.20) 

The  optimal  estimate,  x(t/fj),  is  defined  as  the  sum  of  the  most  recent  nominal 
estimate  and  the  optimal  estimate  of  <5x(f),  denoted  as 

±(t/ti)  =  Xn(t/ti )  +  5±(t/ti).  (2.21) 

Since,  for  the  EKF,  5x(f/fj)  is  zero  over  the  entire  duration  between  measurement 
intervals,  the  best  estimate  of  the  total  state  over  this  interval  would  be  the  solution 
to  the  following 


±{t/ti)  =  f[  x(t/ti),  u(t),  t  ].  (2.22) 

By  calculating  the  partial  derivative  of  f  [  x(f),  u (t),t }  with  respect  to  the  state  vector, 
x,  the  dynamics  partial  derivative  matrix,  F[  t;  x(f/fj)  ],  is  derived  as 
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<9f  [  x(t),  u (t),t 
<9x 


x(t)=x„(t/ti) 


(2.23) 


F[  t ;  x(t/ti)  ] 


It  is  assumed  that  F[  t;  x(f/fj)  ]  is  valid  and  constant  with  respect  to  each  time  inter¬ 
val.  Therefore,  the  transition  matrix  <3>(fj+1,fj)  for  the  stochastic  difference  equation 
can  be  computed  for  each  time  interval  as 


$(ti+1,ti)  =  eF^'  1  At  (2.24) 

5x(ti+i)  =  $(ti+ 1,  ti)  6x(ti)  +  w (ti).  (2.25) 

Similarly,  the  partial  derivative  of  h  [  x(fj),fj  ]  can  be  calculated  from  the 
measurement  equation  with  respect  to  the  state  vector,  x.  The  observation  partial 
derivative  matrix,  H[  U]  x(f“)  ],  is  derived  as 


H[  U;  x(f.  )  ]  = 


<9h  [  x(f),  ti 


dx 


x(t)=x(q  ) 


(2.26) 


Next,  the  estimate  is  propagated  forward  to  the  next  sample  time  tl+\  by  the 
following  estimate  and  covariance  propagation  equations 


x{t/ti)  =  f[  x(t/ti),  u (t),  t  ]  (2.27) 

P(t“+1)  =  <h(fi+i,fi)P(f+)$T(fi+i,fi)  +  (2.28) 

where  Q d(ti)  is  the  discrete  process  noise  intensity.  In  Equation  (2.28),  the  con¬ 
ventional  Kalman  filter  covariance  equation  above  is  used  since  it  is  assumed  that 
F[  f;  x(f/fj)  ]  remains  constant  over  the  time  interval.  Furthermore,  since  this  equa- 
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tion  is  a  relatively  simple  to  calculate  (unlike  the  differential  equation  form  that  will  be 
seen  with  the  other  filters),  the  overall  run  time  of  this  filter  will  be  greatly  decreased. 

Finally,  the  extended  Kalman  filter  measurement  update  equations  incorporate 
the  measurement  by 


K(«0  =  P(i.-)HT[«(;  *(*-)]  (H[4;  x((-)]P(«-)HT[ii;  x(J-)]  +  R(«())_1  (2.29) 

x((,+)  =  x(*r)  +  K(‘<)  (z<  -  h[*(V)>  if])  (2-30) 

P+((i)  =  (I  -  K((i)HT[ij;  x(t-)])  P(*r),  (2.31) 

where  K(£j)  is  the  Kalman  gain  and  I  is  the  identity  matrix.  Now  that  the  new  whole 
value  state  estimate  x.(tf)  is  defined,  it  can  be  used  to  reset  the  new,  most  recent 
nominal  estimate  xn(t/ti+i)  and  is  proceeded  to  the  next  propagation  set  as 

xn(t/ti+1)  =  x(tf)  (2.32) 

x(t/ti+1)  =  xn(t/ti+ 1)  +  6x(t/ti+ 1).  (2.33) 


2.4  Orbital  Mechanics 

Accurately  predicting  an  orbiting  satellite’s  position  at  any  given  instant  is 
vital  if  precise  geolocation  of  ground  targets  is  to  be  achieved.  Satellite  tracking 
can  be  accomplished  by  one  of  two  methods:  an  onboard  navigation  system,  such  as 
INS,  or  an  external  tracking  system  with  an  onboard  receiver,  such  as  GPS.  Orbital 
prediction  is  often  determined  several  hours  (up  to  48  hours)  in  advance  and  is  based 
on  the  known  orbital  parameters  of  the  satellite  [15].  In  order  to  fully  understand  the 
behavior  of  satellites  orbiting  the  Earth,  an  in-depth  explanation  of  orbital  mechanics 
is  warranted,  including  an  understanding  of  Newton’s  laws,  Kepler’s  laws,  Keplerian 
element  transformation,  and  orbit  types. 
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2-4-1  Newton’s  Laws  .  To  understand  how  a  satellite  travels  around  a  planet, 
the  laws  of  motion  derived  by  the  English  scientist  Sir  Isaac  Newton  (1642  -  1727) 
must  be  analyzed  [31].  Based  on  Newton’s  second  law,  that  is,  the  acceleration  of 
the  center  of  a  body  is  proportional  to  the  force  applied  to  that  body,  this  concept 
can  be  used  with  regard  to  the  satellite  and  the  Earth  as  point  mass  bodies  under 
gravitational  attraction.  Although  gravitational  forces  are  not  the  only  forces  acting 
upon  these  bodies,  it  is  by  far  the  largest  contributor  and  is  a  valid  approximation  of 
orbital  position  prediction  [15]. 

According  to  Newton,  the  gravitational  force  of  the  Earth  onto  an  orbiting 
satellite  is 


GMm  r 


(2.34) 


where  G  =  6.67  x  ICE11  N(m/kg)2  is  the  Earth’s  gravitational  constant,  M  =  5.97  x 
1024  kg  is  the  mass  of  the  Earth,  m  is  the  mass  of  the  satellite  in  kg,  r  is  the  distance 
between  the  point  masses  in  km  (equaling  the  radius  of  the  Earth  plus  the  satellite’s 
altitude  above  the  Earth),  and  r  =  rsat  —  r Earth  is  the  position  of  the  satellite  relative 
to  the  Earth,  in  a  Cartesian  coordinate  system,  in  km.  The  motion  of  the  satellite 
can  therefore  be  written  as  a  second  order  differential  equation  as  [19] 


r  + 


G(M  +  m) 

- ^ - r  =  0 


and  can  be  simplified,  since  M  S>  m,  to 


(2.35) 


r  4 — -y- r  =  0  (2.36) 

where  GM  =  3.986  x  108  m3/s3  is  the  Earth’s  standard  gravitational  parameter. 

2-4-2  Kepler’s  Laws  .  It  is  well  understood  that  objects  orbiting  the  Earth 
follow  an  elliptical  pattern.  One  of  first  to  discover  this  was  the  sixteenth-century  Ger- 
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Orbital  Ellipse 

1  St  Law 

Figure  2.7:  Kepler’s  three  laws  of  planetary  motion:  orbital  ellipse,  law  of  equal  areas, 
and  orbits  with  equal  periods  [22], 

man  astronomer  Johannes  Kepler  (1571-1630).  Kepler  devised  the  three  fundamental 
laws  of  planetary  motion  as  follows  [22]: 

•  Planets  follow  an  elliptical  orbit  with  the  sun  at  one  of  its  foci. 

•  A  line  joining  a  planet  and  the  sun  sweeps  out  equal  areas  during  equal  intervals 
of  time. 

•  Two  planets  with  the  same  semi-major  axis  length  have  equal  orbital  periods. 
These  laws  are  summarized  in  Figure  2.7. 

2-4-3  Keplerian  Elements  .  Given  knowledge  of  Kepler’s  laws  of  planetary 
motion  under  idealized  conditions  (i.e. ,  a  perfect  ellipse),  the  motion  of  a  satellite 
can  be  characterized  by  an  elliptical  orbit  in  space  with  the  Earth  as  one  of  its  foci. 
This  orbit  can  be  specifically  identified  in  three  dimensions  by  six  geometrical  orbital 
parameters,  known  as  geometric  Keplerian  elements  [15]  [22],  Five  of  these  elements 
describe  the  size  and  shape  of  the  orbit,  while  the  sixth  element  describes  the  position 
of  the  satellite  at  a  particular  instant  in  time  (or  “epoch”).  The  definitions  of  these 
six  geometric  Keplerian  elements  are  listed  below,  and  are  depicted  in  Figure  2.8: 

•  Semi-major  axis  (a):  The  size  of  the  orbit  in  km,  its  length  is  the  distance 
between  the  geometric  center  of  the  orbital  ellipse  and  the  apoapsis  (point  of 
farthest  approach  to  the  central  body). 


2nd  Law 


P:  Peiiori 

M:  Mojoi  Axis  Lemith 
P3/M3  same  foi  all  planets 

3rd  Law 
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Figure  2.8:  The  six  geometric  Keperian  orbital  elements  in  two  and  three  dimensions, 
respectively.  The  X  and  Y  directions  correspond  to  the  semi-major  and 
semi-minor  axes  of  the  orbit,  respectively.  The  elements  a,  e,  i,  u  and 
v  describe  the  position  of  the  satellite  P  at  a  given  epoch  [22], 

•  Eccentricity  (e ):  The  shape  of  the  orbit.  It  is  a  measure  of  how  much  the  ellipse 
deviates  from  a  perfect  circle  (when  e  =  0). 

•  Inclination  (i):  The  angle  between  the  orbital  plane  to  central  body’s  equator 
in  rad. 

•  Right  Ascension  of  the  Ascending  Node  (RAAN  or  Rt)\  The  rotation  of  orbit’s 
reference  plane  with  respect  to  ascending  node  (the  point  on  the  satellite’s  orbit 
where  it  crosses  the  equatorial  plane)  in  rad. 

•  Argument  of  Perigee  (u):  The  angle  from  the  ascending  node  to  perigee  (the 
point  at  which  the  satellite  is  closest  to  the  center  of  the  Earth)  in  rad. 

•  True  Anomaly  (v):  The  location  of  the  satellite  at  a  given  epoch,  this  is  the 
angle  between  the  ascending  node  and  the  satellite  position  in  the  orbital  plane 
in  rad. 

Given  the  Keplerian  elements,  the  satellite’s  orbital  period  in  sec,  Torb,  and  orbital 
apogee  in  m,  aporb,  (the  point  in  orbit  farthest  from  the  Earth)  can  be  deduced  as 
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Torb  =  2,71 


(2.37) 


a-Porb  —  (1  +  e)a  —  Re,  (2.38) 

where  GM  is  the  standard  gravitational  parameter  (see  Section  2.4.1)  and  Re  is  the 
radius  of  the  Earth  in  m. 

2-4-4  Satellite  Position  and  Velocity  Transformation  .  Although  it  is  more 
common  to  see  the  Keplcrian  elements  described  in  Section  2.4.3  as  geometrical  char¬ 
acterizations  of  the  orbit  itself,  these  elements  also  can  be  represented  as  the  six 
satellite  position  and  velocity  elements  in  the  Cartesian  coordinate  system  [15].  The 
defined  orbital  coordinate  system  is  equivalent  to  the  inertial  frame  (i-frame)  in  which 
the  origin  is  located  at  the  elliptical  foci  located  at  the  center  of  the  Earth.  The  x 
and  y  axes  of  this  reference  frame  correspond  to  the  semi-major  and  semi-minor  axes 
of  the  orbit,  respectively,  as  seen  previously  in  Figure  2.8. 

In  order  to  generate  the  three-dimensional  satellite  position  and  velocity  vectors, 
r  and  v,  in  the  i-frame,  the  six  geometric  Keplerian  elements  a,  e,  i,  hi,  u  and  v  must 
undergo  a  reference  frame  transformation  [15].  The  magnitude  of  the  orbital  radius, 
|| r j| ,  or  the  distance  between  Earth’s  center  and  the  satellite,  is  defined  as 

a(l  —  e2) 

1  +  e  cos  v 

The  satellite  position  and  velocity  vectors  in  the  two-dimensional  orbital  plane  de¬ 
scribed  in  Figure  2.8  (denoted  here  as  the  pqw- frame),  are  calculated,  respectively, 
as 


(2.39) 
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(2.40) 


and 


Transforming  from  the  pqw- frame  to  the  i- frame,  the  direction  cosine  matrix  Cpqv, 
(see  Section  2.2.3)  is  defined  as 


cos  12  cos  u  —  sin  12  sin  u>  cos  i  —  cos  12  sin  u)  —  sin  12  cos  u)  cos  i  sin  Q  sin  i 

Gpqw  =  sin  12  cos  oj  —  cos  12  sin  u>  cos  i  —  sin  12  sin  u>  +  cos  12  cos  u>  cos  i  —  cos  12  sin  i 

sin  uj  sin  i  cos  u  sin  i  cos  i 

(2.42) 

Finally,  the  satellite  position  and  velocity  vectors  in  the  i- frame  are  calculated,  re¬ 
spectively,  to  be 


24 


and 


s’ 'ii 

^ pqw 


rpqw 


(2.43) 


v  = 


=  C 


pqw 


pqw 


(2.44) 


2.4.5  Orbit  Types  .  Now  that  the  elliptical  orbit  of  the  satellite  can  be  fully 
characterized  as  the  six  satellite  position  and  velocity  parameters,  a  unique  elliptical 
orbit  can  be  designed.  Considerations  for  the  satellite’s  orbit  are  largely  based  upon 
the  desired  orbital  period  and  apogee  (see  Section  2.4.3).  For  comparison,  three 
popular  orbital  classifications  are  listed  below  and  depicted  in  Figure  2.9  [22,31]: 


•  Low  Earth  Orbit  (LEO):  The  orbit  closest  to  the  Earth’s  surface.  It  lies  just 
beyond  the  thermosphere  (or  outer  atmosphere).  LEO  orbits  typically  have  a 
period  of  approximately  90  min  and  an  apogee  between  450  and  600  km  above 
the  Earth’s  surface.  High-resolution  imaging  satellites  and  the  International 
Space  Station  (ISS)  are  located  in  this  orbit. 

•  Geosynchronous  Orbit  (GEO):  A  high  Earth  orbit  (HEO),  a  geosynchronous 
orbit  has  a  period  equal  to  that  of  Earth’s  (approximately  24  hrs).  Its  apogee 
is  approximately  36,000  km  above  the  Earth’s  surface.  A  GEO  orbit  directly 
above  the  Earth’s  equator  is  known  as  a  geostationary  orbit.  Communication, 


25 


Figure  2.9:  Three  common  satellite  orbits:  low  Earth  orbit  (LEO),  Geosynchronous 
orbit  (GEO),  and  Molniya  orbit  (MOL),  the  last  two  of  which  are  defined 
as  high  Earth  orbits  (HEOs)  [31]. 

early  warning  and  nuclear  detection  satellites  are  located  in  geosynchronous  and 
geostationary  orbits. 

•  Molniya  Orbit  (MOL):  The  Molniya  orbit,  a  HEO  orbit  named  after  the  Sovi¬ 
et/Russian  satellites  of  the  same  name,  has  a  unique  12  hr  orbit  at  an  approxi¬ 
mate  63  deg  inclination.  With  an  apogee  between  26,000  and  38,000  km  above 
the  Earth’s  surface,  a  distinct  advantage  of  satellites  in  MOL  orbits  is  their 
ability  to  track  locations  in  the  northern  hemisphere  for  an  extended  period  of 
time,  with  less  ground  interference  at  high  look  angles  than  would  be  true  with 
lower  Earth  orbits.  Communication  and  intelligence  satellites  use  this  orbit. 

2.5  Optical  Modeling  and  Spatial  Coherence 

In  this  section,  the  basic  physical  properties  of  an  optical  sensor  model  are 
presented,  and  the  concept  of  spatial  coherency  is  reviewed. 

2.5.1  Optical  Modeling  .  An  optical  sensor  is  designed  to  measure  the  inten¬ 
sity  of  light  entering  the  device  through  an  aperture  and  converts  it  into  an  electrical 
signal  which  can  be  read  by  an  observer  [25,28].  These  sensors  typically  make  use 
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of  either  a  charge-coupled  device  (CCD)  or  complimentary  metaloxide  semiconductor 
(CMOS)  active-pixel  sensor  to  create  a  two-dimensional  image  as  a  function  of  light 
intensity  [9]. 

For  the  purposes  of  this  research,  the  world  is  defined  as  a  collection  of  real 
objects  of  interest  that  illuminate  the  world  and  interact  with  the  other  physical 
objects  through  various  types  of  reflection.  In  radiometry,  irradiance  is  defined  as  the 
amount  of  light  intensity  with  respect  to  unit  area  that  falls  on  an  object’s  surface  [2], 
The  physical  irradiance  pattern  enters  the  aperture  of  the  optical  sensor  (defined  as 
the  scene)  and  is  projected  onto  the  image  plane.  This  process  is  represented  as 
a  continuous  array  of  nonnegative  real  numbers.  For  simplicity,  object  surfaces  are 
assumed  to  be  Lambertian,  meaning  the  brightness  of  the  surface  to  an  observer  is 
the  same  regardless  of  the  observer’s  angle  of  view. 

A  digital  optical  imaging  sensor  consists  of  an  aperture,  lens,  detector  array,  and 
sampling  array.  A  simple  imaging  system  model  is  shown  in  Figure  2.10.  The  lens 
images  the  scene  on  the  detector  array.  The  light  pattern  focused  on  the  detector  array 
is  defined  as  the  image.  The  detector  array  converts  the  light  energy  into  a  voltage 
or  a  charge  which  is  converted  to  a  digital  value  by  the  sampling  array  (e.g.,  an  8-bit 
digital  value  within  the  set  [0-255],  where  255  represents  the  highest  intensity). 

2.5.2  Spatial  Coherence.  Spatial  coherence  is  defined  as  the  property  of 
waves  to  maintain  definite  phase  in  space.  Within  the  topic  of  imaging,  an  under¬ 
standing  of  coherence  is  warranted  as  it  describes  the  correlation  properties  of  light 
waves  [9] .  The  constructive  addition  and  destructive  subtraction  of  light  waves,  known 
as  optical  interference,  can  affect  the  resulting  image  clarity,  particularly  in  the  pres¬ 
ence  of  atmospheric  turbulence. 

First,  the  physical  properties  of  light  waves  must  be  characterized.  In  mathe¬ 
matical  terms,  the  equation  for  a  light  wave  is  defined  as  a  complex  held 

Ug{u,  v )  =  A{u ,  v)e~i  ^u’v)  (2.45) 
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Figure  2.10:  A  digital  optical  imaging  system,  consisting  of  an  aperture,  lens,  detec¬ 
tor  array  and  sampling  array  [28].  Objects  in  the  world  are  collected 
into  the  aperture  (a  scene)  as  an  irradiance  pattern  and  focused  onto 
the  detector,  producing  an  image  represented  digitally  as  an  array  of 
nonnegative  numbers. 


where  A(u,  v )  is  the  amplitude,  v )  is  the  phase  and  i  is  the  imaginary  number. 

Coherent  sources  have  correlated  (non-random)  wave  phase  while  incoherent  sources 
have  uncorrelated  (random)  phase  [8]. 

The  simulation  of  coherent  and  incoherent  sources  through  a  lens  can  be  achieved 
by  the  means  of  convolution  [8].  This  convolution,  in  which  the  object  forms  a  spread 
function  region  at  the  image  plane,  results  in  a  slightly  blurred  image.  This  method 
is  depicted  in  Figure  2.11. 

It  can  be  shown  that  coherent  imaging  is  linear  in  amplitude  while  incoherent 
imaging  is  linear  in  intensity,  as 


Ii(u,  v)  —  \h(u,  v)  ®  Ug(u,  v)\2  (Coherent  Imaging) 
Ii(u,  v)  —  \h(u,  v)\2  ®  Ig(u,  v)  (Incoherent  Imaging), 


(2.46) 

(2.47) 


Figure  2.11:  Convolution  through  a  lens.  An  object  at  point  P  forms  a  spread  function 
of  region  Q,  resulting  in  a  slightly  blurred  image. 

where  h(u,v )  is  the  amplitude  spread  function,  \h(u,v)\2  is  the  point  spread  function, 
Ig(u,v)  is  the  object  intensity,  /*(«,  v)  is  the  image  intensity,  and  “(8)”  is  the  convo¬ 
lution  operator  [8].  For  coherent  light,  image  intensity  is  the  squared  convolution  of 
the  object  wave  function  and  the  amplitude  spread  function.  For  incoherent  light,  im¬ 
age  intensity  is  the  convolution  of  the  object  intensity  and  the  point  spread  function 
(which  is  the  amplitude  spread  function  squared).  The  resulting  effects  of  these  cases 
are  illustrated  with  a  simple  rectangular  object  pattern  in  Figure  2.12  [21],  In  the 
case  of  coherent  imaging,  the  resulting  image  intensity  has  a  blurred  “waffling  effect” , 
a  phenomenon  known  as  Fresnel  ringing,  which  is  the  result  of  optical  diffraction  [21], 
In  the  case  of  incoherent  imaging,  the  image  intensity  distributes  evenly,  producing 
images  with  sharper  features  than  those  produced  by  coherent  imaging. 

In  the  case  of  satellite  imaging,  the  sun  is  a  predominantly  incoherent  source, 
meaning  the  majority  of  the  light  rays  from  the  sun  have  random  phase.  This  results 
in  nearly  uniform  phase  distribution  on  the  image  plane  as  light  from  the  sun  bends 
around  the  edges  of  the  object  (or  objects)  in  the  object  plane.  Therefore,  image 
blurring  as  the  result  of  Fresnel  ringing  is  not  of  concern. 
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Figure  2.12:  Coherent  and  incoherent  imaging  with  a  rectangular  pattern.  The  object 
(a)  is  convolved  with  the  amplitude  (or  point  spread)  function  (b)  to 
produce  either  a  coherent  image  (c)  or  incoherent  image  (d).  Note  the 
blurred  “waffle  effect”  in  (c),  the  result  of  optical  diffraction  [21], 
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2.6  Imaging  Through  Turbulence 

Ground  images  are  continuously  captured  from  space  by  orbiting  observation 
satellite  systems.  When  these  images  are  collected  from  a  satellite  several  kilometers 
above  the  Earth’s  surface,  it  is  likely  that  image  resolution  and  image  displacement 
will  be  affected  by  the  atmosphere.  In  order  to  understand  the  turbulent  effects  on  an 
image  collected  through  the  atmosphere,  a  background  of  turbulence  theory  must  be 
presented.  Following  this  background,  the  concepts  and  derivations  of  the  turbulence 
strength  parameter,  C%  is  defined.  The  atmospheric  coherence  width,  r0,  among  other 
important  turbulence  parameters,  will  then  be  discussed  in  light  of  the  problem  at 
hand.  Unless  stated  otherwise,  all  equations  in  Section  2.6  are  derived  in  Ref.  [1], 

2.6.1  Turbulence  Background  .  The  Earth’s  atmosphere  is  comprised  of 
gases,  chemicals  and  water  vapor,  all  contributing  to  refractive  index  fluctuations, 
causing  light  waves  to  bend  and  scatter  (or  attenuate)  unpredictably.  The  refractive 
index  of  the  Earth’s  atmosphere  is  near  unity,  however,  as  light  propagates  through 
this  medium,  the  optical  waves  become  randomly  distorted  and  the  resulting  image 
resolution  becomes  limited  [21]. 

Atmospheric  turbulence  is  the  result  of  stochastic  variations  in  temperature 
and  velocity  within  the  Earth’s  atmosphere  [21].  This  turbulence  is  caused  by  a 
combination  of  solar  heating  the  surface  and  atmosphere,  convection  (causing  hot  air 
to  rise  and  cold  air  to  fall)  and  diffusion  (mixing  areas  of  high  concentration  and 
areas  low  concentration).  An  illustration  of  the  atmospheric  layers  is  depicted  in 
Figure  2.13.  The  largest  concentration  of  turbulence  is  within  the  first  20  to  24  km 
above  the  surface  [1].  It  is  assumed  in  this  analysis  that  weather  conditions  are  clear 
with  no  obstruction  due  to  clouds,  rain  or  fog. 

2. 6.2  Turbulence  Strength  .  The  index  of  refraction,  n,  is  one  of  the  most  sig¬ 
nificant  parameters  with  respect  to  light  wave  propagation  through  the  atmosphere  [1] . 
Fluctuations  in  atmospheric  refractive  index  are  related  to  corresponding  fluctuations 
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Figure  2.13:  The  Earth’s  atmospheric  layers.  Nearly  100%  of  the  atmospheric  mass 
is  contained  within  the  troposphere  and  stratosphere,  the  majority  of 
which  is  contained  solely  within  the  troposphere. 

in  temperature,  T  (measured  in  K),  and  pressure,  P  (measured  in  millibars),  and  can 
be  written  accordingly  as 


n(R)  =  l  +  79xlO-6^,  (2.48) 

where  R,  assigned  as  any  point  in  space,  is  bounded  by  the  inertial  subrange,  or  the 
range  of  unstable  air  masses  (eddies)  defined  by  the  inner  scale  bound  lQ  and  outer 
scale  bound  La  as  [lQ  -C  R  -C  L0\.  For  this  research,  the  Kolmogorov  turbulence 
spectrum  will  be  assumed;  therefore,  the  inertial  subrange  is  unbounded. 

The  statistical  description  of  the  random  turbulence-induced  fluctuations  in  the 
atmosphere’s  refractive  index  can  be  expressed  as  a  structure  function  with  respect 
to  R.  Assuming  statistically  homogeneous  and  isotropic  turbulence  ( R  =  |Ri  —  R2I2), 
the  structure  function,  Dn(R),  is  expressed  as 


Dn{R) 


ic2nlo4'3  R2  0  <R</0 
[C2R2/ 3  Z0<R<L0, 


(2.49) 


where  C2  is  the  refractive  index  structure  parameter.  C2  is  a  measure  of  the  strength  of 
fluctuations  in  the  refractive  index,  and  can  be  interpreted  as  a  measure  of  the  strength 
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of  atmospheric  turbulence.  C2n  is  quantified  in  terms  of  pressure,  P,  temperature,  T 
and  the  temperature  structure  constant,  Cf.,  by 


C2n  =  (j9  x  KT6^  C2.  (2.50) 

The  associated  power  spectral  density  for  refractive  index  fluctuations,  over  the 
inertial  subrange  defined  by  [l/f„<K€l/4  is  defined  by 

$„(«)  =  0.0033C'2kT11/3.  (2.51) 

Cl  is  known  to  vary  as  a  function  of  height  above  ground,  the  strongest  occurring 
during  the  daytime  near  the  ground  (on  an  order  of  10-13m-2/3  or  higher).  In  order 
to  express  the  turbulent  strength  as  a  function  of  height,  models  for  both  C2  and 
atmospheric  winds  must  be  assigned.  The  Button  wind  model  is  commonly  used  to 
describe  the  atmospheric  winds  as 


V (z)  =  vg  +  usz  cos  (z  +  vT  exp 


z  COS  (z  —  Ht 


2  Al1/2 


sin2  0  +  cos2  0  cos2  ( z\ 


(2.52) 


where  z  is  the  propagation  slant  path  (in  m)  given  by  the  equation  z  —  h  sec  where 
h  is  the  height  above  ground  in  m  and  (  is  the  zenith  angle  in  rad.  The  zenith  angle 
angle  assignment  is  depicted  in  Figure  2.14.  Furthermore,  Vq  =  5  m/s  is  assigned 
as  the  ground  wind  speed,  u>s  =  0  deg/s  as  the  turbulent  slew  rate,  vt  =  30  m/s  as 
the  tropopause  wind  speed,  Ht  =  12.5  km  as  the  average  altitude  of  the  tropopause, 
Lt  =  4.8  km  as  the  average  thickness  of  tropopause  and  0  =  0  deg  as  the  direction 
of  wind  speed  [1]. 

Using  the  Bufton  wind  model  in  Equation  (2.52),  the  root  mean  square  (RMS) 
wind  speed  due  to  turbulence  is  calculated  to  be 
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Figure  2.14:  Zenith  and  nadir  viewing  angles,  (z  and  (n,  respectively.  The  relation¬ 
ship  of  slant  path,  z,  with  respect  to  altitude,  h,  is  z  =  hsec(  [21]. 
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(2.53) 


where  the  integration  limits  in  Equation  (2.53)  encompass  the  troposphere  defined  in 
Figure  2.13. 

Finally,  the  Hufnagel- Valley  (H-V)  5-7  model,  commonly  used  to  describe  C2, 
is  defined  as  [1] 


C2n(h)  =  0.00594  ( 
+2.7  x  IQ”16 


^)2(10-5ft)'»exp 

eXP(i^)+"leXP 


(2.54) 


where  h  =  zcosC,,  A  =  1.7  x  10-14  m _2/3,  and  vrms  =  21  m/s  from  Equation  (2.53). 
Note  that  A  is  the  nominal  value  of  C2  at  the  ground  (or  C^(0)).  Using  the  H-V 
C2  model,  a  comparison  of  the  turbulent  strength,  C2  as  a  function  of  propagation 
distance,  z,  and  zenith  angle,  (,  can  be  seen  in  Figure  2.15. 
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Figure  2.15:  Turbulence  strength,  C%,  as  a  function  of  propagation  distance,  z,  and 
zenith  angle,  (Z-  As  (z  increases,  the  relative  height  above  ground  de¬ 
creases  to  where  turbulence  is  greatest,  contributing  to  higher  for 
longer  durations. 
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2.6.3  Atmospheric  Coherence  Width  .  The  atmospheric  coherence  width, 
ro,  also  called  Fried’s  parameter,  is  the  measure  of  the  spatial  coherence  of  light  at 
the  receiver.  A  large  ro  implies  that  light  is  coherent  across  large  distances,  indicating 
good  imaging  system  performance  [1,21].  Since  turbulent  strength,  C%,  is  nonlinear  in 
nature,  the  direction  of  propagation  can  greatly  affect  the  resulting  coherence  width. 

Consider  two  distinct  paths  of  a  propagating  point  source,  producing  a  spherical 
wave  at  the  source  that  becomes  a  planar  wave  over  a  large  propagation  distance 
(such  as  the  case  between  the  Earth’s  surface  and  an  orbiting  satellite).  In  the  first 
case,  an  uplink  path  is  defined  in  which  the  atmosphere  is  farthest  from  the  receiver 
(i.e.,  an  imaging  system  in  orbit).  In  the  second  case,  a  downlink  path  is  defined  in 
which  the  turbulent  atmosphere  is  nearest  the  receiver  (i.e.,  for  an  imaging  system 
on  the  ground).  The  resulting  coherence  width  of  the  uplink  path,  onto  the  satellite 
receiver,  is  many  times  larger  than  the  satellite  itself.  This  large  coherence  width 
is  due  to  the  sun’s  incoherent  light,  reflected  from  the  ground,  bending  appreciably 
early  in  the  propagation  path  (when  the  wave’s  spatial  extent  is  small)  and  remaining 
essentially  unchanged  as  it  propagates  toward  the  satellite.  The  opposite  is  true  for 
the  relatively  small  coherence  width  of  the  downlink  path,  in  which  the  wave’s  spatial 
extent  is  large  when  entering  the  turbulent  atmosphere  and  remains  turbulent  as  it 
propagates  toward  the  ground  [6,7].  An  illustration  of  these  two  cases  can  be  seen  in 
Figure  2.16. 

For  a  satellite  imaging  system,  in  which  light  is  propagated  from  the  ground 
(the  transmitter)  to  the  satellite  (the  receiver),  the  uplink  path  is  assigned.  For 
coherence  width  computation,  it  is  necessary  to  identify  the  statistical  path  moments 
of  the  optical  system.  The  uplink  propagation  first  and  second  path  moments  are, 
respectively, 
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l^2u 


[H  c&h)  [e  +  e{]5/3dft 

J  h0 

r  ci(h)(i~  tf 3  dh, 

Jho 


(2.55) 

(2.56) 


where  h  is  the  the  elevation  in  km,  ho  =  0  is  the  elevation  of  the  ground  (at  sea  level), 
H  is  the  elevation  of  the  satellite  in  km,  and  £  =  (h  —  ho)/(H  —  ho)  is  the  normalized 
distance  variable.  The  output  (receiver)  beam  curvature  parameter,  0,  equals  one  for 
a  spherical  wave.  Therefore,  0  =  1  —  0  =  0.  The  resulting  coherence  width  for  the 
uplink  path,  r'o„  .  is  defined  as 


[0.42  sec(Gv)*;2  (tii u  +  0.622/^A11/6)]  3/5 

(2.57) 

[0.42 sec((N)k'2idlu]~3/5 . 

(2.58) 

Since  the  output  (receiver)  Fresnel  ratio  parameter,  A,  equals  zero  for  a  spherical  wave 
case,  Equation  (2.57)  is  simplified  to  Equation  (2.58).  Note  that,  for  the  uplink  case, 
the  viewing  angle,  £jv,  is  with  respect  to  nadir  at  the  satellite,  pointing  downward.  As 
expected,  r0  for  the  uplink  case  will  be  much  larger  than  for  the  downlink  case.  This 
result  implies  good  imaging  system  performance  for  an  observation  satellite  collecting 
images  of  the  ground. 

2.6.4  Image  Jittering  .  An  on-axis  object  may  appear  to  wander  within 
a  satellite  image  due  to  atmospheric  turbulence.  This  frequent  shifting  (or  displace¬ 
ment)  of  the  object,  otherwise  known  as  “image  jittering”  is  associated  with  the  angle 
of  arrival  fluctuations  of  an  optical  wave  onto  the  receiver  aperture  [1,32].  These  angle 
of  arrival  fluctuations,  <  (3%  >,  are  defined  primarily  by  turbulence  at  high  altitudes, 
and  are  related  to  the  coherence  width  of  the  propagation  path  (see  Section  2.6.3). 


37 


Figure  2.16:  Uplink  and  downlink  path  scenarios  for  a  propagating  point  source 
(spherical  wave).  Note  the  effects  of  free-space  diffraction  in  the  uplink 
case,  providing  minimal  atmospheric  turbulence  effects  and  ultimately  a 
large  coherence  width  at  the  receiver  (the  satellite)  [21]. 


The  RMS  angle  of  arrival  for  an  uplink  path,  in  rad,  simplified  for  a  spherical 
wave,  is  calculated  as 


<  ft  >=  2.91/iiusec(Gv)(JDG)-1/3,  (2.59) 

where  Dq  is  the  aperture  diameter  of  the  receiver  in  m.  The  resulting  image  displace¬ 
ment  standard  deviation  (or  “image  jitter”),  in  m,  is  calculated  to  be 

CTimag  =  flens\J<  PI  >,  (2.60) 

where  fiens  is  the  focal  length  of  the  lens  in  m.  The  variance  of  the  angle  of  arrival 
fluctuations  in  rad,  a2,  which  attributes  to  reduction  in  image  resolution,  is  calculated 
as 


aa  =  6-88(r0)  5/3sec(C N)(DG)  1/3k  2,  (2.61) 
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where  r0  is  the  uplink  coherence  width  in  m  and  k  is  the  wave  number  in  rad/m. 
In  the  uplink  case,  crirnag  and  are  very  small,  since  r0  is  very  large.  Assuming  an 
off-nadir  angle  of  less  than  25  deg  for  LEO  orbits,  and  an  angle  of  less  than  5  deg  for 
HEO  orbits,  atmospheric  turbulence  dne  to  image  jittering  can  generally  be  ignored 
with  respect  to  an  observation  satellite. 

2.6.5  Atmospheric  Refraction  .  In  free  space,  light  waves  propagate  in 
straight  lines  due  to  the  fact  that  their  dielectric  permittivity,  eo,  and  magnetic  per¬ 
meability,  po,  are  constant  in  space  and  time  relative  to  the  speed  of  wave  propagation 
c  as  defined  by  the  following  formula  [8] 


c  =  Gu0e0)-1/2.  (2.62) 

However,  the  dielectric  permittivity  of  the  atmosphere,  e,  is  greater  than  that  of  free 
space  (e  >  e0).  Therefore,  light  waves  propagate  at  a  speed  v  that  less  than  c  and  in 
doing  so  deviate  from  straight  propagation  paths,  resulting  in  refraction  (or  bending) 
of  the  beam  [5].  This  results  in  a  non-uniform  atmospheric  refraction,  n,  that  is  greater 
than  unity  with  respect  to  the  Cartesian  coordinate  x—,  y —  and  z— directions,  or  [8] 

n(x,y,z )  =  -  >  1.  (2.63) 

Therefore,  in  order  to  determine  the  true  path  of  the  light  wave,  atmospheric  refraction 
must  be  considered. 

According  to  Ref.  [1],  the  refractive  index  of  the  atmosphere  can  be  closely 
approximated  as  a  function  of  optical  wavelength,  A  (in  m),  pressure,  P  (in  mbars), 
and  temperature,  T  (in  K).  Assuming  the  non-uniform  lower  atmospheric  density  is 
limited  to  within  20  km  of  the  surface,  P  and  T  are  analyzed  as  a  function  of  altitude, 
z  (in  km).  From  [5],  pressure  in  the  lower  atmosphere  is  calculated  to  be 
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(2.64) 


P(z)  =  ||  -  (2.26  x  10~5z  -  l)5'32xl°3||,  0  <  z  <  20  km. 

Likewise,  temperature  in  the  lower  atmosphere  is  computed  as 

{-6.81  x  10~3(^  -  1)  +  293,  0  <  *  <  11  km 

(2.65) 

218,  11  <  2  <  20  km. 

Finally,  the  atmosphere’s  refractive  index,  n,  can  be  written  in  terms  of  optical  wave¬ 
length,  A  (in  m),  pressure,  P  (in  millibars),  and  temperature,  T  (in  K),  as  [1] 


n(z) 


1  +  77.6  x  10“6(1  +  7.52  x  10“3A"2) 


P(z) 
T(z)  ’ 


0  <  z  <  20  km. 


(2.66) 


Figure  2.17  depicts  the  plot  of  pressure,  temperature  and  refractive  index,  with  respect 
to  altitude. 

Next,  Snell’s  law  [9]  is  used  to  calculate  the  ground-to-satellite  horizontal  dis¬ 
placement  due  to  light  refraction.  Snell’s  law  is  defined  as 


ri\  sin  6*i  =  n2  sin  d2  =  n3  sin  93,  (2.67) 

where  n i ,  n2  and  n3  are  the  refractive  indices  of  the  troposphere,  tropopause  and 
free  space,  respectively.  This  relationship  can  be  seen  in  Figure  2.18.  Given  a  known 
off- zenith  departure  angle  6 1,  the  sequential  angles  9-2  and  03  can  also  be  computed. 
Using  geometry,  the  horizontal  displacement  Axn  can  be  calculated  as 


Axn 


Xw  XyjQ 

(zw i  tan  9W i  +  zw2  tan  0w2  +  zw3  tan  9w3)  -  zw0  tan  9, 
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(2.68) 

(2.69) 
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Figure  2.17:  Pressure,  temperature  and  refractive  index,  with  respect  to  altitude,  at 
a  wavelength  of  675  nm.  Between  11-20  km,  temperature  remains  steady 
through  the  tropopause.  Above  an  altitude  of  20  km,  it  is  assumed  that 
turbulence  is  negligible,  and  therefore  is  approximately  where  free  space 
begins  (n  ~  1). 
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Figure  2.18:  The  refraction  of  a  light  wave  through  the  tropopause  (n2)  to  free  space 
(n3).  Since  the  index  of  refraction  decreases  to  approximately  1  at  alti¬ 
tudes  above  the  tropopause  (n3  <  n2),  light  waves  refract  to  a  greater 
extent,  resulting  in  a  larger  departure  angle  in  free  space  (03  >  d-2 ) 

where  xw  is  the  total  horizontal  displacement  from  the  surface  to  the  satellite  including 
atmospheric  refraction  and  xwo  is  the  horizontal  displacement  from  the  surface  to  the 
satellite  without  atmospheric  refraction.  Assuming  an  off-nadir  angle  of  less  than  25 
deg  for  LEO  orbits,  and  an  angle  of  less  than  5  deg  for  higher  orbits,  the  effects  of 
atmospheric  refraction  is  found  to  be  minimal. 

2.7  Image- Aided  Navigation 

In  previous  research  [25,27],  it  has  been  demonstrated  that  the  coupling  of 
imaging  with  inertial  sensors  has  provided  a  navigation  improvement  of  at  least  two 
orders  of  magnitude  over  inertial  systems  without  the  aid  of  optical  devices.  As 
such,  knowledge  of  the  navigation  state  (i.e.,  position,  velocity,  and  attitude)  of  a 
space  vehicle  could  also  be  improved  by  using  image-aided  navigation,  provided  that 
optical  measurements  from  an  imaging  sensor  pointed  toward  the  ground  are  available. 
This  section  will  summarize  the  background  behind  image-aided  navigation,  review 
a  means  of  matching  a  satellite  image  to  predefined  image  template,  and  introduce 
the  concept  of  georeferencing  in  order  to  identify  a  targets  location  on  a  relative 
coordinate  system. 
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2.7.1  General  Background  .  Although  advances  in  the  field  of  image-aided 
navigation  have  been  made,  the  level  of  accuracy  in  such  a  system  is  critically  de¬ 
termined  by  the  alignment  and  calibration  of  the  imaging  sensor  [23].  Previous  ap¬ 
proaches,  including  mechanical  techniques  and  field-calibrated  estimation  based  tech¬ 
niques  [11],  have  had  limited  success,  requiring  dedicated  equipment  unsuitable  for 
field  work  and  being  subject  to  intermittent  manufacturer  errors. 

Other  more  recent  approaches  use  real-time  estimators  and  use  the  field  of  visible 
stars  to  provide  the  reference  for  the  optical  system  [29].  No  operator  involvement 
or  external  equipment  is  required  in  this  stellar  observation  approach.  Additionally, 
star  observation  accounts  for  time-varying  errors  prominent  in  inertial  sensors  (see 
Section  2.2.4)  and  has  the  advantage  of  operating  in  real  time.  An  example  of  such 
a  spacecraft  navigation  system  includes  the  high  accuracy  star  tracker  (HAST)  [14]. 
However,  stellar  observations  require  visibility  of  the  sky  and  star  tracking  algorithms 
must  be  sensitive  enough  to  resolve  the  location  of  celestial  objects.  Additionally,  if 
an  onboard  imaging  system  has  the  combined  role  of  star  tracking  as  well  as  Earth 
observation,  the  ability  to  track  both  star  and  ground  targets  may  be  limited. 

2.7.2  Image  Matching  .  Image  matching  is  used  in  many  applications,  in¬ 
cluding  object  recognition,  stereo  matching,  and  feature  tracking,  as  a  means  of  iden¬ 
tifying  a  feature  or  area  common  among  a  series  of  pixelated  images  [3].  Image-aided 
navigation  can  benefit  from  image  matching  algorithms  by  incorporating  matched  im¬ 
age  data  between  subsequent  images  to  determine  how  a  stationary  target  “wanders” 
between  those  images. 

Image  matching  applications  used  for  feature  tracking  commonly  use  the  sum- 
of-squared- difference  (SSD)  measure  to  determine  the  best  match  between  subsequent 
images  [16].  The  SSD  method  operates  directly  on  the  pixelated  irradiance  pattern 
of  the  image  (see  Section  2.5.1)  and  measures  the  correspondence  between  an  image 
and  a  template  (in  which  a  match  of  the  original  image  exists).  Knowledge  of  the 
maximum  position  uncertainty  of  the  feature  of  interest  within  the  template  allows 
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the  image  matching  to  be  performed  efficiently  with  greatest  likelihood  of  determining 
a  unique  image  match. 

The  SSD  can  be  computed  as  [3] 

n  n 

SSD(u,v )  =  EE  [f(x,  y)  -  t(x  -  u  +  dx,y  -  v  +  dy)\ 2  (2.70) 

x=0  y= 0 

where  /  is  the  satellite  image  and  t  is  the  template  in  which  the  image  is  to  be  matched 
within.  The  summation  is  over  the  Cartesian  coordinates  (x,y),  corresponding  to  the 
northern  and  eastern  directions  in  the  n-frame,  respectively.  A  template  with  northern 
and  eastern  boundary  limit  parameters,  u  and  v,  form  a  “window”  with  respect  to 
the  origin  (the  origin  is  assumed  to  be  in  the  center  of  the  image).  Disparities  dx 
and  dy  are  the  vertical  and  horizontal  differences  between  images  in  the  n-frame, 
respectively,  and  are  sometimes  used  to  determine  depth  or  distance  of  an  object 
within  the  image  [17].  Low  SSD  values  represent  a  good  match  between  the  image 
and  template,  indicating  where  the  difference  between  the  image  and  template  is 
minimal.  In  Equation  (2.70),  the  differencing  between  the  image  and  template  is 
squared  to  ensure  that  SSD  results  are  always  non-negative. 

Although  the  sum-of-squared-difference  has  the  advantage  of  being  relatively 
simplistic  and  easy  to  implement,  weaknesses  of  SSD  also  exist.  This  method  is  sen¬ 
sitive  to  outliers  and  is  not  accustomed  to  template  variations,  such  as  those  that 
occur  at  occluding  boundaries  in  the  image  [16].  Since  SSD  is  restricted  to  over¬ 
lapping  the  entire  image  over  the  template,  this  algorithm  can  be  easily  “fooled” 
by  repeated  patterns  throughout  the  template,  resulting  in  increasing  probability  of 
erroneous  matching.  This  weakness  is  especially  inherent  to  indoor  navigation,  in 
which  repeatable  features  (i.e.,  flooring  and  ceiling  tiles)  are  common.  Fortunately, 
this  problem  is  not  as  prevalent  with  outdoor  navigation  and  satellite  imagery,  where 
feature  repeatability,  particularly  at  higher  resolutions,  is  rare. 
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Figure  2.19:  Georeferenced  control  points  with  unique  features  are  used  to  georegis¬ 
ter  a  target  in  a  satellite  image  to  a  reference  map  (image  provided  by 
Manifold  Systems). 

2. 7.3  Georeferencing  .  Georeferencing  is  the  process  of  identifying  a  target’s 
location  on  a  coordinate  system  relative  to  the  Earth.  The  process  of  georeferencing 
a  target  onto  a  reference  map  is  depicted  in  Figure  2.19,  where  pre-assigned  feature- 
rich  control  points  within  the  image  are  used  to  “georegister”  (or  adjust)  a  target  to 
a  geographic  location  based  upon  a  known  reference  map.  This  process  is  known  as 
georegistration.  The  target’s  georeferenced  location  (measured  in  e-frame  latitude/- 
longitude  coordinates)  is  referred  to  as  the  target’s  geolocation. 

A  georeferenced  image  can  be  used  to  measure  the  coordinates  of  a  target  of 
interest  within  an  image.  However,  inaccurate  geolocation  can  limit  this  capability, 
for  example,  if  an  insufficient  number  of  control  points  are  available  to  accurately 
determine  the  target’s  location.  Large  satellite  attitude  errors  characterized  by  a 
low-performance  inertial  sensors  can  also  limit  geolocation  performance. 

Tightly  coupled  image-aided  inertial  navigation  systems  have  been  designed  to 
extract  navigation  information  of  the  satellite  by  automatically  detecting  and  tracking 
stationary  optical  features  of  opportunity  in  the  environment,  thereby  vastly  reducing 
the  vehicle’s  attitude  errors  [26].  One  significant  advantage  of  this  navigation  system 
is  that  it  can  operate  in  areas  where  GPS  is  either  denied  or  unavailable. 
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2.8  Stochastic  Constraints 


Using  an  imaging  sensor  to  determine  the  navigation  states  from  optical  mea¬ 
surements  depends  upon  tracking  the  target  location  through  a  series  of  images.  Tar¬ 
get  recognition  through  image  interpretation,  however,  can  be  very  inefficient  and  time 
consuming  if  the  number  false  matches  are  not  limited  and  the  feature  correspondence 
search  is  not  constrained.  Therefore,  the  stochastic  projection  method  developed  in 
Refs.  [28,30],  which  constrains  this  correspondence  search  by  area  by  incorporating  a 
priori  knowledge  of  the  satellite  navigation  states,  allows  the  user  to  accurately  and 
optimally  predict  the  pixel  location  and  uncertainty  of  a  target  feature  in  a  series  of 
images. 

As  mentioned  in  Ref.  [30],  this  stochastic  projection  method  uses  many  of  the 
assumptions  of  the  Kalman  filter  in  Section  2.3.1.  The  landmark  of  interest  is  assumed 
to  be  stationary  (or  very  slowly  moving)  with  respect  to  the  surface  of  the  Earth. 
Additionally,  the  camera  is  assumed  to  be  rigidly  mounted  to  the  vehicle  with  a  known 
alignment  and  calibration.  Finally,  it  is  assumed  that  the  terrain  is  flat  (constant 
elevation).  All  equations  in  Section  2.8  are  derived  in  Ref.  [30]. 

Given  the  navigation  state  at  time  tt,  x(U),  described  in  Equation  (2.4),  the 
landmark  position  corresponding  to  a  pixel  location,  y(U),  is  a  non-linear  function 
of  x,  given  by 


y(U)  =  g  [x(U)].  (2.71) 

The  calculated  landmark  position,  y(U),  is  modeled  as  a  perturbation  about  the  true 
position  as 


y(U)  =  y(U)  +  £y(*i),  (2-72) 

where  y(U)  is  a  function  of  the  calculated  navigation  state,  x(fj).  x(fj)  is  also  modeled 
as  a  perturbation  about  truth,  and  is  of  the  form 
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St(u)  =  x(tj)  +  <5x(tj). 


(2.73) 


Applying  perturbation  techniques  to  the  landmark  position  function,  the  land¬ 
mark  error,  Sy(ti),  can  be  expressed  as  a  linear  function  of  the  navigation  state  errors 


Sy  (U)  =  G  yx(ti)6x(ti), 


(2.74) 


where  the  influence  coefficient,  G yx(ti),  is  defined  as 


Gyx(ti) 


dg  [x(U)', 


d  x(tj 


(2.75) 


x(ti)=x(tt) 


The  landmark  error  covariance,  P yy(ti),  and  cross-correlation.  P xy(ti),  are  defined  as 


p  yyiti)  =  E[5  y(ti)hyr(f?:)]  (2.76) 

Pxy(U)  =  E[6x(ti)5yT(ti)}.  (2.77) 

Substituting  Equation  (2.74)  into  Equation  (2.76)  yields 

P  yy(ti)  =  Gyx(ti)Pxx(ti)GyX(ti).  (2.78) 


The  cross  correlation  matrices  are  calculated  in  a  similar  manner  as 


P  xy(ti)  =  P„(ti)G  IM  (2.79) 

P  yxiti)  =  Gyx(ti)Pxx(ti).  (2.80) 
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Combining  the  navigation  state  with  the  landmark  state  through  augmentation, 
the  initial  combined  state  vector,  x*(i0),  and  initial  combined  covariance  matrix, 
P*(tG),  are  defined,  respectively,  as 


x*(tc)  = 


x(t0 

y(G 


and 


P*(«o)  = 


Pa :xi^o)  Pa;y(to) 

Pyx  (tQ)  Pyy(t0) 


Likewise,  the  combined  psuedonoise  matrix,  G *(U),  is  expressed  as 


G  x(ti)  0 
0  G  y(ti 


and  the  combined  process  noise  intensity,  Q *(ij),  is  expressed  as 


Q  *(U)  = 


Q  X(U)  o 

o  Q  y(ti 


G y(ti)  dehnes  the  landmark  error  dynamics  as  a  random  walk  by 


(2.81) 


(2.82) 


(2.83) 


(2.84) 


Sy  (ti)  =  Gy(ti)wy(ti), 


(2.85) 


and  w y(ti)  is  a  zero-mean,  white  Gaussian  noise  process  with  covariance  kernel 
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E[wy(ti)Wy(ti  +  t)]  =  Q  y(ti)S(T). 


(2.86) 


For  propagating  the  combined  navigation  and  landmark  estimate,  the  EKF 
nonlinear  stochastic  differential  equation  is  expressed  as 


x*(t/U)  =  f  [  x*(t/ti),  u*(t),  t  }  +  G*(ti)w (ti)>  (2.87) 

in  which  the  dynamics  of  the  landmark  state  within  x*(t/ti)  are  zero.  Therefore,  the 
respective  partial  derivative  dynamics  matrices  for  the  navigation  and  landmark  state 
are  derived  as 


F*[  t;  ii(t/ti)  }  = 


df  [  x,ux(t),t 


dx 


(2.88) 


X=Xn(t/tj) 


Fy[  t;  y (t/U)  }  = 


df  [  y,u y(t),t' 


dy 


(2.89) 


y=y  n(*Ai)=° 


and  the  combined  partial  derivative  dynamics  matrix  is  written  as 


F*[  t;  rit/U)  }  = 


Fx[t;  x(t/U ) 


0  Fy[  t\  y(t/U)  ] 


(2.90) 


The  pixel  location  measurements  for  the  EKF  are  modeled  as  a  known  nonlinear 
pixel  projection  function  of  the  combined  state  plus  linearly  additive  measurement 
noise  as 


z (ti)  =  h  [  x*(ti),ti  ]  +v(ti), 


(2.91) 
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such  that  the  pixel  location  measurements  are  a  function  of  the  combined  navigation 
and  landmark  state.  The  respective  partial  derivative  pixel  projection  matrices  are 
then  derived  as 


tt  r  .  I  A  <9h  [  X,ti 


<9x 


X=x(ti  ) 


(2.92) 


Hu/[  U]  y {ti )  ]  = 


a  cdi  [  y,  ^ 


<9y 


(2.93) 


y=y  (*»  ) 


and  the  combined  partial  derivative  pixel  projection  matrix  is  written  as 


H*[  U-  x*(tr)  ]  = 


H2X[  til  x(ti  )  ]  0 


0  Hzy[  tj;  y(ti  )  ; 


(2.94) 


Finally,  the  the  covariance  of  the  pixel  location  errors  can  be  written  as  the 
EKF  residual  covariance  computed  within  the  Kalman  gain  (Equation  (2.30))  as 


P zz(ti)  =  H*[  ti;  x*(tr)  ]P*(4r )H*T[  tg  x*(fr)  ]  +  R^).  (2.95) 

Therefore,  given  the  pixel  coordinates  of  a  stationary  ground  landmark  at  time 
tj,  the  predicted  pixel  coordinates  of  the  same  landmark  at  time  tj+i  can  be  described 
by  the  pixel  location  error  covariance,  Pzz(tj),  as  a  function  of  the  combined  navigation 
and  landmark  state  covariance,  P*(tj),  the  pixel  projection  matrix,  H*[  t*;  x*(t~)  ], 
and  the  measurement  noise  intensity,  R (ti). 


50 


III.  Method 


In  this  chapter,  the  modeling  and  methodology  of  a  satellite-based  image-aided 
navigation  system  will  be  covered,  using  concepts  covered  in  the  previous  chap¬ 
ter.  First,  two  satellite  trajectories,  one  low-Earth  orbit  and  one  high-Earth  orbit, 
will  be  produced.  Next,  modeling  of  the  satellite  system  will  be  developed,  includ¬ 
ing  the  vehicle’s  imaging  system  parameters  as  well  as  identifying  system-level  noise 
parameters.  The  ground  image  model  will  be  constructed,  using  the  defined  satellite 
trajectory  and  noise  parameters,  as  well  as  existing  satellite  imagery.  Image  match¬ 
ing  and  georeferencing  techniques  will  be  implemented  to  predict  the  landmark  state. 
Finally,  an  extended  Kalman  filter  model  will  be  presented. 

3.1  Orbit  Modeling 

One  of  the  first  steps  in  developing  an  image-aided  navigation  system  for  an 
orbiting  satellite  is  the  understanding  of  the  satellite’s  trajectory  around  the  Earth. 
This  requires  knowledge  of  the  satellite’s  position  and  velocity,  both  of  which  are 
computed  from  this  trajectory. 

In  order  to  prove  orbital  independence  of  this  satellite-based  system,  two  orbit 
types,  one  LEO  and  one  HEO  (specifically,  a  MOL),  will  be  produced.  Descriptions 
of  these  orbits  can  be  found  in  Section  2.4.5.  Using  geometric  Keplerian  orbital 
elements  defined  in  Section  2.4.3,  both  orbits  are  characterized.  The  basis  of  this 
element  assignment  may  be  dependent  upon  orbit  requirements  for  the  system,  such 
as  a  specified  orbital  period  or  orbital  apogee  (see  Equations  (2.37)  and  (2.38)).  It 
should  be  noted  that  these  orbit  declarations  are  assumed  to  be  nominal  and  are  not 
corrupted  by  error.  The  parameters  assigned  for  the  low  Earth  and  Molniya  orbits  of 
this  satellite  are  represented  in  Table  3.1. 

3.2  Satellite  System  Modeling 

Provided  the  assigned  Keplerian  orbit  elements  mentioned  in  Section  3.1,  the 
initial  conditions  for  all  navigation  and  landmark  states  of  this  satellite  navigation 
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Table  3.1:  Assigned  nominal  Keplerian  elements  for  a  typical  LEO  and  MOL.  Using 
these  elements,  the  orbital  period  and  apogee  were  calculated. 


Keplerian  Element 

LEO 

MOL 

a  (km) 

6,760 

26,600 

e 

0.0323 

0.704 

i  (rad) 

1.71 

1.10 

fl  (rad) 

5.97 

3.60 

c o  (rad) 

5.46 

4.92 

v  (rad) 

3.12 

1.36 

Torb  (hrs) 

1.54 

12.0 

aporb  (km) 

598 

38,900 

system  can  be  defined.  These  nominal  parameters  will  later  provide  a  basis  for  the 
true  navigation  and  landmark  states  (where  error  corruption  will  be  introduced). 

3.2.1  Nominal  System  State  Modeling  .  First,  the  initial  position  and 
velocity  navigation  states  for  the  nominal  satellite  system  dynamics  are  defined  in 
the  three-dimensional  Cartesian  coordinate  system.  Using  Equations  (2.39)  through 
(2.44),  presented  in  Section  2.4.4,  the  position  vector  in  the  i-frame,  rx,  ry ,  and  rz  (in 
units  of  km)  and  velocity  vector  in  the  i-frame,  rx ,  ry ,  and  rz  (in  units  of  km/min), 
are  computed  from  the  geometric  Keplerian  elements  presented  in  Table  3.1. 

Next,  the  initial  attitude  states  of  the  satellite  in  the  n-frame,  cj)  and  9,  that  being 
the  pitch  and  roll  angles  of  the  satellite  imaging  system  measured  in  rad,  are  assumed 
to  nominally  be  zero.  In  other  words,  it  is  assumed  that  initially  the  satellite’s  imaging 
sensor  is  pointed  directly  at  the  target  of  interest  on  the  ground.  For  simplicity,  it 
is  assumed  that  the  x-  and  y-axes  of  the  satellite  image  (in  the  n-frame)  are  always 
aligned  with  respect  to  the  true  latitude  and  longitude  axes  of  the  Earth,  respectively; 
therefore,  image  rotation  is  not  of  concern  and  the  yaw  angle  of  the  satellite  need  not 
be  estimated. 
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Table  3.2:  Initial  nominal  LEO  and  MOL  navigation  and  landmark  states  at  apogee, 
with  respect  to  their  respective  i  and  n-frames.  The  initial  navigation 
states,  rXnom0 ,  rVnorn0,  rZnom0 ,  vXnom0,  vVnom0,  and  vZnom0,  are  derived  using 
Keplerian  orbital  element  transformation.  Assuming  the  satellite’s  imaging 
sensor  is  pointed  directly  at  the  target,  the  initial  attitude  states,  4>nOm0  and 
6nom0,  and  the  initial  landmark  states,  tn Xnom0  and  tnyrlorn0,  are  considered 
to  be  zero. 


Initial  Conditions 

LEO 

MOL 

rxnom0  (km) 

-4,710 

17,100 

fynom  0  (km) 

780 

-13,900 

rznom0  (km) 

5,080 

39,600 

vXnom0  (km/min) 

-300 

81.2 

Vynomo  (km/min) 

139 

50.0 

Vznomo  (km/min) 

-299 

-17.3 

(ftnomO  (l  sd) 

0 

0 

Onomo  (rad) 

0 

0 

tnxnomO  (pixels) 

0 

0 

tnynomO  (pixels) 

0 

0 

Finally,  the  initial  n- frame  landmark  states  for  the  system  are  defined  as  the 
northern  and  eastern  target  position  errors  on  the  ground,  tnx  and  tny.  Since  this  state 
measurement  is  based  upon  the  resolution  of  the  satellite  image,  it  is  calculated  in 
units  of  pixels  (see  Section  2.5.1  for  details  diffraction-based  imaging).  Similar  to  the 
attitude  states,  the  nominal  landmark  states  are  initially  assumed  to  be  zero,  again 
implying  the  satellite’s  imaging  sensor  is  pointed  directly  at  the  target  of  interest  on 
the  ground.  For  the  derived  low  Earth  and  Molniya  orbits,  the  satellite  navigation 
and  landmark  states  at  apogee  are  represented  in  Table  3.2. 

3.2.2  Satellite  System  Dynamics  .  The  navigation  and  landmark  state  dy¬ 
namics  equations  can  now  be  presented.  Rewriting  Equation  (2.36)  from  Section  2.4.1, 
the  satellite’s  acceleration  vector  in  the  i- frame,  r,  can  be  computed  as 
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r  = 


(3.1) 


GM 


3 


r  +  wf  (t), 


where  GM  is  the  Earth’s  standard  gravitational  parameter  (approximately  398,600 
km3 /min2),  ||  r  ||  is  the  magnitude  of  the  satellite’s  position  vector  and  w, (t)  is 
uncorrelated,  zero-mean,  white,  Gaussian  noise  source  with 


E[wr(t)wr(t)T(t  +  r)]  =  qr(t)6(r).  (3.2) 

The  randomized  displacement  of  the  satellite  image  in  the  n-frame  as  a  result 
of  variations  in  the  attitude  of  the  imaging  sensor  (assume  to  be  rigidly  fixed  to  the 
satellite  body)  can  be  described  as  a  first-order  time-correlated  drift.  This  drift  can 
be  characterized  as  a  first-order  Gauss-Markov  process  [12],  in  which  the  satellite’s 
rate  of  pitch  and  rate  of  yaw,  respectively,  can  written  as 


— </>  +  w^(t) 

T<t> 

(3.3) 

—  6+Wg(t), 

Tg 

(3.4) 

where  and  Tg  are  the  first-order  Gauss  Markov  time  constants  of  the  roll  and  pitch 
angle  in  arcmin,  respectively,  and  w ^(t)  and  w g(t)  are  uncorrelated,  zero- mean,  white, 
Gaussian  noise  sources  with 


E[w^(t)w^t)T(t  +  r)]  =  q^(t)S(r) 

(3.5) 

E[wd(t)wd(t)T{t  +  t)]  =  qg(t)6(r). 

(3.6) 

As  explained  in  Section  2.8,  the  landmark  state  vector  will  be  used  to  accurately 
predict  the  pixel  location  of  a  target.  As  such,  tnx  and  tny  will  be  defined  as  a 
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Figure  3.1:  Generated  nominal  satellite 
trajectory  in  a  LEO. 


Figure  3.2:  Generated  nominal  satellite 
trajectory  in  a  MOL. 


stochastic  projection  of  the  navigation  states,  and  will  therefore  have  no  dynamics. 
In  other  words, 


tn  =  0  +  w  in(t),  (3.7) 

where  wt^  (f )  is  uncorrelated,  zero-mean,  white,  Gaussian  noise  source  with 

^KWwt„WT(i  +  r)]  =  9tn(f)W-  (3.8) 

Using  the  initial  conditions  found  in  Table  3.2  and  the  defined  satellite  system 
dynamics  expressed  in  Equations  (3.1)  through  (3.7),  the  nominal  low-Earth  and 
Molniya  orbits  are  generated.  The  resulting  path  of  the  satellite’s  position  in  each 
orbit  is  depicted  in  Figures  3.1  and  3.2. 
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Table  3.3:  Assigned  imaging  parameters  for  an  observation  satellite,  for  both  LEO 
or  MOL  orbits.  A  high  apogee  in  the  MOL  orbit  results  in  higher  image 
resolution,  PX,  a  smaller  maximum  off-nadir  angle,  r]Nmaxi  and  smaller 
view  angle  (the  maximum  angle  at  which  the  image  can  be  viewed).  Since 
similar  imaging  systems  will  be  used  for  either  orbit,  the  operating  wave¬ 
length,  A,  lens  diameter,  Diens  and  lens  focal  length,  fiens,  will  not  vary. 


Imaging  Parameters 

LEO 

MOL 

PX  (krn/pix) 

0.1 

1.0 

(Nmax  (deg) 

25 

5.0 

view  angle  (deg) 

0.4 

0.2 

A  (nm) 

675 

675 

Diens  (jft) 

0.60 

0.60 

fiens  (ffl) 

10 

10 

It  is  assumed  that  the  optimal  operation  time  of  the  satellite  is  when  it  reaches 
its  orbital  apogee,  since  it  is  at  that  location  in  orbit  where  the  satellite  travels  at 
its  slowest  rate  (as  dictated  by  Kepler’s  2nd  law  in  Section  2.4.2)  and  therefore  the 
satellite  is  able  to  capture  as  many  images  of  the  target  of  interest  as  possible  within 
the  designated  off-nadir  angle  limits. 


3.2.3  Image  System  Modeling  .  Based  on  an  appropriate  imaging  system 
suitable  for  either  LEO  or  MOL  orbits  [4,18,20],  the  observation  satellite’s  assigned 
imaging  parameters  are  represented  in  Table  3.3. 

It  will  be  assumed  that  the  designed  LEO  orbit  with  an  apogee  of  598  km  will 
have  an  image  resolution,  PX,  of  less  than  1  km/pixel  (perhaps  on  the  order  of  0.1 
km/pixel),  whereas  the  MOL  orbit  with  an  apogee  of  38,900  km  will  likely  have  a 
PX  of  no  better  than  1  km/pixel.  Likewise,  the  satellite’s  off-nadir  angle,  r]Nmax,  and 
maximum  angle  at  which  the  image  can  be  viewed,  view  angle,  will  likely  be  much 
larger  for  a  satellite  in  a  LEO  orbit  versus  that  in  a  MOL  orbit,  is  the  satellite’s 
location  is  orders  of  magnitude  closer  to  the  Earth’s  surface  in  the  lower  orbit.  For 
simplicity,  it  will  be  assumed  that  all  other  imaging  system  parameters  would  be 
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sufficient  for  either  orbit.  Therefore,  its  operating  wavelength,  A,  lens  diameter,  Diens 
and  lens  focal  length,  fiens,  will  remain  unchanged. 

3.3  Noise  Modeling 

It  is  unreasonable  to  assume  that  the  nominal  image-aided  satellite  navigation 
system  described  in  Section  3.2  would  not  be  corrupted  by  real-world  errors.  Such 
errors  could  include  measurement  noise  clue  to  sensor  misalignment  (Section  2.2),  tra¬ 
jectory  noise  due  to  atmospheric  drag  (Section  2.4),  image  sensor  noise  clue  to  vehicle 
vibration  (Section  2.5),  or  image  displacement  (or  fluctuation)  clue  to  atmospheric 
turbulence  (Section  2.6). 

3.3.1  Turbulence  Noise  .  Errors  clue  to  atmospheric  turbulence  are  identified 
in  order  to  determine  if  these  errors  are  significant  enough  to  model  in  the  satellite 
system.  As  discussed  in  Section  2.6,  it  is  already  assumed  that  these  errors  will  be 
minimal  in  even  the  worst  case  scenarios  (i.e. ,  when  the  off- nadir  angle  (n  =  C »„„)• 
For  completeness,  however,  these  calculations  are  warranted. 

From  Sections  2.6.4  and  2.6.5,  it  was  shown  that  target  location  errors  due  to 
turbulence  may  exist  due  to  image  jitter  and  horizontal  light  refraction  displacement. 
Using  Table  3.3,  and  Equations  (2.60)  and  (2.69),  the  resulting  image  jitter,  <Jimg,  and 
refraction  displacement,  Axn,  are  summarized  in  Table  3.4. 

As  expected,  at  minimum  off- nadir  angles  ((n  =  0),  turbulence  error  due  to 
image  jitter  and  refraction  displacement  are  insignificant  in  either  orbit.  At  maximum 
off-nadir  angles  (Cat  =  C Nmax ),  &img  is  again  insignificant,  and  Axn  is  approximately 
0.5  pixels  (with  respect  to  image  resolution)  for  either  orbit.  It  can  be  deduced  that 
for  uplink  propagation  at  relatively  high  altitudes  above  the  lower  atmosphere  (i.e., 
the  case  of  an  imaging  sensor  onboard  an  orbiting  satellite),  image  jitter  is  not  found 
to  be  of  concern,  and  image  displacement  due  to  refraction  is  only  of  slight  concern 
when  off-nadir  angles  approach  their  maximum  limit.  As  such,  for  completeness,  very 
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Table  3.4:  Image  position  error  due  to  turbulence  in  LEO  and  MOL  orbits.  At  mini¬ 
mum  off-nadir  angles  (Cat  =  0),  image  jitter,  aimg,  and  refraction  displace¬ 
ment,  Axn,  are  practically  non-existent  for  either  orbit.  At  maximum 
off-nadir  angles  (Ov  =  C Nmax),  &img  is  again  very  small,  and  Axn  is  only  a 
half  a  pixel  (with  respect  to  image  resolution)  in  either  case. 


Error  due  to  Turbulence 

LEO 

Molniya 

Tmg  at  (N  =  0  (pixels) 

4.8  xlO-12 

1.0  xlO"34 

a*mg  at  c n  =  C Nmax  (pixels) 

5.0  xlO-12 

1.1  xlO-34 

Axn  at  Cat  =  0  (pixels) 

0 

0 

Axn  at  C n  =  C Nmax  (pixels) 

0.54 

0.56 

slight  additive  error  will  be  introduced  in  the  attitude  and  landmark  states,  as  listed 
in  the  following  section. 


3.3.2  Image  Sensor,  Trajectory  and  Measurement  Noise  .  As  mentioned 
earlier  in  Section  3.3,  space  vehicles,  as  well  as  their  onboard  hardware,  are  subject  to 
errors  as  a  result  of  real-world  noise  sources,  including  vehicle  vibration,  atmospheric 
drag  and  measurement  miscalibration,  just  to  name  a  few.  In  order  to  model  these 
errors  in  a  recursive  estimator  (such  as  the  extended  Kalman  filter  derived  in  Sec¬ 
tion  2.3.2),  these  noise  sources  are  all  assumed  to  be  zero- mean,  Gaussian  and  white. 
For  the  purpose  of  modeling,  these  assumptions  are  assumed  to  be  reasonable.  With 
this  in  mind,  careful  consideration  must  be  given  as  to  the  numerical  assignment  of 
these  noise  sources. 

The  assigned  initial  statistical  and  process  noise  parameters  for  the  designed 
LEO  and  MOL  orbits  are  summarized  in  Tables  3.5  and  3.6,  respectively.  Given 
faster  satellite  trajectories,  larger  off-nadir  angles  and  higher  atmospheric  effects  in 
a  LEO  orbit,  it  is  reasonable  to  assume  that  all  initial  uncertainties  be  an  order  of 
magnitude  smaller  and  all  process  noise  strengths  be  at  least  an  order  of  magnitude 
larger  than  they  would  be  in  a  MOL  orbit.  An  exception  would  be  the  time  constants, 
T0  and  To,  since  the  rate  of  image  sensor  drift  is  independent  to  the  altitude  of  the 


Tabic  3.5:  Assigned  initial  statistical  parameters  for  a  satellite  in  LEO  and  MOL 
orbits.  Faster  satellite  trajectories,  larger  off-nadir  angles  and  higher  at¬ 
mospheric  effects  in  the  LEO  orbit  justify  these  parameters  to  be  an  order 
of  magnitude  smaller  than  they  would  be  in  the  MOL  orbit,  with  the  ex¬ 
ception  of  the  independent  time  constants,  ly,  and  tq. 


Statistical  Parameters 

LEO 

MOL 

^VxOf  ®ry0)  @rz 0  (km) 

0.5 

0.05 

Vvyo,  crVz0  (km/min) 

0.05 

0.005 

0>o,  O0o  (rac0 

5  xlO-4 

5  xlO-5 

re  (min) 

10 

10 

^tnxmeas i  & tnymeas  (pix) 

0.5 

0.05 

Table  3.6:  Assigned  process  noise  parameters  for  a  satellite  in  LEO  and  MOL  orbits. 

Faster  satellite  trajectories,  larger  off-nadir  angles  and  higher  atmospheric 
effects  in  the  LEO  orbit  justify  these  parameters  to  be  at  least  an  order  of 
magnitude  larger  than  they  would  be  in  the  MOL  orbit. 


Noise  Parameters 

LEO 

MOL 

Qrx,  qry,  (-lfz  (km2 /min2) 

5  xlO"4 

5  xlO-6 

%  (min"1) 

5  xlO-8 

5  xlO”10 

Qiny  (pix2/min) 

10 

0.1 

satellite.  Note  that  all  initial  statistical  parameters  are  relatively  small  since  it  is 
assumed  that  a  navigation  system  calibration  had  just  recently  occurred. 

It  should  be  noted  that  the  satellite  roll  and  pitch  rate  noise  strengths,  q ^ 
and  qg,  are  both  first-order  Gauss-Markov  processes  as  described  in  Section  3.2.2. 
Accordingly,  these  values  are  calculated  as 


%  =  2^  (3-9) 

m  =  2=£.  (3,10) 
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Figure  3.3:  Generated  true  satellite  attitude  error,  as  first-order  Gauss-Markov  pro¬ 
cesses  [12].  In  this  particular  example,  an  imaging  system  in  a  MOL  orbit 
is  modeled  over  a  24  hour  period.  As  a  first-order  Gauss-Markov  process, 
the  mean  diverges  toward  zero  at  a  rate  relative  to  the  time  constants 
and  re. 


The  results  of  this  first-order  Gauss-Markov  process,  that  being  the  true  satellite 
attitude  states,  (fttme  and  9true,  are  plotted  in  Figure  3.3.  The  nominal  roll  and  pitch 
are  assumed  to  be  zero  for  all  time  (with  no  optical  drift). 

Furthermore,  the  combined  process  noise  intensity  matrix,  Q*,  including  both 
navigation  and  landmark  states,  is 
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Likewise,  the  measurement  noise  intensity  matrix,  R,  in  pixels,  is  described  as 


(3.12) 


The  pixel  location  measurement  in  the  n-frame,  a  function  of  both  the  navigation  and 
landmark  states,  is  calculated  as 


(3.13) 


and  where  the  slant  range  in  the  e-frame,  dQ ,  from  the  satellite  to  the  ground  target 
coordinates  (xoe,  yoe,  zoe),  is  computed  as 
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%oe)  H-  (j*y  Voe)  z  %oe)  • 


(3.14) 


3.4  Ground  Error  Modeling 

Given  the  assigned  satellite  trajectory  and  noise  parameters,  a  ground  image 
model  can  be  constructed  using  existing  satellite  imagery.  A  three-hundred  square 
mile  image  of  the  greater  Cincinnati,  Ohio  area,  provided  by  the  National  Aeronautics 
and  Space  Administration  (NASA)  Landsat  satellite  [10],  is  assigned  to  represent 
the  area  of  the  Earth  where  the  satellite  is  tracking  the  target  of  interest,  ft  is 
assumed  that  the  satellite  is  at  or  near  its  apogee  so  that  the  target  tracking  duration 
is  maximized.  In  order  to  accurately  represent  the  capabilities  of  a  satellite  image 
sensor  in  either  LEO  or  Molniya  orbits,  the  resolution  of  the  image,  PX,  must  be 
considered.  For  simplicity,  identical  images  of  the  area  are  used  for  either  orbit,  one  in 
which  PX  =  0.1  km/pixel  (or  10  pixels/km)  to  represent  the  satellite  in  a  low  Earth 
orbit,  and  one  where  PX  =  1  km/pixel  to  represent  the  satellite  in  a  Molniya  orbit. 
A  sample  of  the  image  template  (the  area  in  which  the  subsequent  satellite  images 
will  be  captured  within),  at  a  resolution  of  1  km/pixel,  is  depicted  in  Figure  3.4. 

In  order  to  generate  the  truth  data,  for  the  landmark  state,  tnxtrue  and  tnytrue ,  a 
full  understanding  on  the  landmark  error  is  necessary.  Rationally,  target  error  should 
be  a  function  of  the  trajectory  of  the  satellite  (as  the  satellite  moves  across  the  sky 
above  the  target)  as  well  as  a  function  of  the  time-correlated  attitude  errors  affecting 
the  image  sensor  (as  described  in  Section  3.3.2).  This  combination  will  be  assumed 
in  order  to  produce  the  desired  truth  data  for  the  landmark  state. 

The  northern  and  eastern  true  target  position  errors  of  the  satellite  in  the  71- 
frame,  tnXorb  and  tnyorb ,  solely  as  a  function  of  orbital  drift,  are  defined  as 
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Figure  3.4:  A  300  square  mile  Landsat  image  of  the  greater  Cincinnati,  Ohio  area  [10]. 

For  the  LEO  case  (where  the  satellite  is  close  to  the  Earth),  an  image 
resolution,  PX,  of  0.1  km/pixel  (10  pixels/km)  is  used,  whereas  for  the 
Molniya  orbit  case  (where  the  satellite  is  very  far  from  the  Earth  at  it 
apogee),  a  PX  of  1  km/pixel  is  used. 
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(3.15) 

(3.16) 


tnxorb  di, 

tnyorb  hf 


A  l  at 


PX 

Alon  . 
«-p^cos(Zatr 


where  Re  is  the  radius  of  the  Earth  in  the  e-frame  (approximately  6,378  km),  and  the 
angular  latitude  and  longitude  separation  in  the  e-frame  (as  a  result  of  the  difference 
between  the  nominal  and  true  satellite  trajectories),  A/of  and  A  Ion,  are  defined  as 


A  lat 

latnom  lat  true 

(3.17) 

Alon 

lOTlnorn  lOTltrue • 

(3.18) 

Since  it  is  assumed  that  the  terrain  is  flat,  knowledge  of  the  true  altitude  error  is  not 
necessary  and  therefore  is  not  computed.  Next,  the  true  target  location  error  as  a 
function  of  only  the  attitude  error  first  order  Gauss-Markov  processes,  tnxFOGM  and 
tnyFOGM,  described  in  Section  3.3.2,  are  calculated  as 


tnxFOGM  ~  d0(j)pOGM  (3.19) 

tnyFOGM  =  da9  FOGM  •  (3.20) 

The  slant  range,  dQ,  from  the  satellite  position  to  the  target  aimpoint,  is  calculated 
as 


do  \J  ( r'x  -t  oe)  T  ( l  y  Hoe )  ~  T  (Pz  /-oe)  i 


(3.21) 


and  the  first-order  Gauss-Markov  attitude  errors  are  calculated  as 
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4>fogm{}) 

^fogm(*) 


(3.22) 

(3.23) 


_  e  A*/r*  (j)FOGM  (i  —  1) 

_  e  At /re  $FOGM(i  —  1), 

where  %  is  the  current  sample,  i  —  1  is  the  previous  sample,  and  At  is  the  time  between 
samples. 

Finally,  the  combined  true  target  location  error  in  the  n-frame,  tnxtrue  and  tnytrue, 
as  a  function  of  both  orbital  drift  and  first  order  Gauss-Markov  processes,  are  defined 
as 


tnx  truth  tnxorb  A  tnxFOGM 

tny  truth  tnyorb  tnyFoGM 


(3.24) 

(3.25) 


An  example  of  the  resulting  landmark  state  truth  data,  mapped  onto  a  target 
of  interest  (East  Fork  Lake,  25  miles  east  of  downtown  Cincinnati),  is  depicted  in 
Figure  3.5. 


3.5  Image  Matching  Development 

Image  matching  techniques  discussed  in  Section  2.7.2  can  be  implemented  as  a 
means  of  tracking  a  target  with  unique  features  among  a  series  of  pixelated  images.  By 
incorporating  matched  image  data  between  subsequent  images,  it  can  be  determine 
how  a  stationary  target  appears  to  “drift”  between  those  images,  therefore  providing 
indirect  knowledge  of  how  the  true  landmark  error  is  generated. 

Using  the  sum-of-squared-difference  approach  defined  in  Equation  (2.70),  a 
satellite  image,  representing  what  the  onboard  image  sensor  “sees”  based  upon  the 
satellites  view  angle,  is  sampled  across  each  pixelated  row  and  column  with  respect 
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Figure  3.5:  Generated  ground  error  in  the  LEO  case,  as  a  function  of  orbital  drift  and 
first-order  Gauss-Markov  processes.  In  (a),  the  target  of  interest  (East 
Fork  Lake)  is  tracked  in  a  generated  image  with  a  swath  width  of  500 
pixels  (50  km),  as  dictated  by  the  satellite’s  view  angle.  The  blue  asterisk 
denoting  the  target  in  (b)  is  within  3  sigma  of  the  generated  ground  error 
uncertainty  (represented  here  as  a  dashed  circle). 
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Figure  3.6:  Matching  a  satellite  image  containing  a  target  of  interest  (b)  within  an 
image  template  (a),  using  the  SSD  approach.  In  (a),  the  overlaying  image 
is  sampled  across  each  row  and  column,  and  the  SSD  for  each  sampling  is 
computed  in  order  to  determine  the  lowest  SSD  intensity  (best  match). 

to  an  image  template  (that  being  the  300  square  mile  image  presented  in  Figure  3.4). 
This  concept  is  represented  in  Figure  3.6,  where  the  target  of  interest  (the  lake  tracked 
in  Figure  3.5)  is  matched  onto  the  300  square  mile  map  of  the  greater  Cincinnati,  Ohio 
area. 

The  SSD  process  can  be  represented  as  a  three  dimensional  contour  plot,  in 
which  lower  SSD  values  indicating  where  the  difference  between  the  image  and  tem¬ 
plate  is  minimal,  representing  the  best  match  between  the  two  images.  This  is  illus¬ 
trated  in  Figure  3.7. 

3.6  Georeferencing  Development 

Although  the  image  matching  process  described  in  Section  3.5  is  an  efficient 
means  of  determining  how  a  target  of  interest  wanders  between  subsequent  satellite 
images,  it  has  no  knowledge  of  any  existing  initial  measurement  error  generated  prior 
to  the  image  matching  process.  In  other  words,  when  the  first  image  is  captured,  the 
image  sensor  assumes  that  no  initial  target  location  error  exists,  and  unknowingly 
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Figure  3.7:  A  3-D  contour  plot  of  the  resulting  SSD  calculations  for  each  row  and 
column  sample.  The  “best  match”  is  determined  to  be  where  the  SSD 
values  reach  a  local  minimum,  indicating  where  the  difference  between 
the  satellite  image  and  image  template  is  smallest. 

matches  the  future  captured  images  relative  to  the  first  image.  Whether  or  not  the 
image  does  in  fact  contain  the  target  in  never  truly  determined. 

Georeferencing  provides  a  solution  to  this  dilemma.  Knowledge  of  the  initial 
target  location  error  by  means  of  automated  georeferencing  (see  Section  2.7.2),  the 
target’s  true  location  in  the  rt- frame  can  be  determined  relative  to  a  reference  map  of 
the  Earth.  Essentially,  a  georeferenced  image  can  be  used  to  accurately  predict  the 
initial  true  pixel  location  of  the  target.  The  benefits  of  geolocation  are  apparent  in 
Figure  3.8,  where  a  constant  “bias”  between  the  SSD  matched  landmark  state  and 
the  truth  data  can  be  seen.  The  remaining  bias  between  the  matched  data  and  truth 
data  is  simply  calculated  as 


t 


nbias 


~  t 


nmatch 


- 1 


ntruth  • 


(3.26) 
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Figure  3.8:  Plotting  of  the  best  matched  northern  and  eastern  target  location  errors 
against  the  generated  truth  data.  The  matched  data  accurately  follows 
the  truth  data,  however,  a  constant  error  remains  because  the  initial  true 
pixel  location  is  unknown.  This  “bias“  between  the  matched  data  and 
truth  data  can  be  corrected  using  georeferencing  techniques. 


Note  that  tnbias  is  only  an  approximation  of  this  bias,  and  not  a  true  representation, 
since  its  calculation  is  limited  by  the  available  image  resolution  (rounded  to  the  nearest 
pixel). 

Through  combined  image  matching  and  georeferencing  techniques,  an  accurate 
target  location  error  prediction  is  generated  that  closely  matches  the  true  target  lo¬ 
cation  error 


t npred 
t npred 


—  t nm.atch  t nbias 

—  Untruth- 


(3.27) 

(3.28) 
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3.7  Extended  Kalman  Filter  Development 

Provided  the  truth  data  for  the  navigation  and  landmark  state  in  Section  3.2, 
the  statistical  and  noise  parameters  generated  in  Section  3.3,  and  the  truth  (and 
predicted)  data  derived  in  Sections  3.4  through  3.6,  the  extended  Kalman  filter  defined 
in  Section  2.3.2  can  now  be  built. 

The  initial  navigation  and  landmark  state  mean,  Xq,  and  state  covariance,  Pq, 
both  derived  in  Equations  (2.81)  and  (2.82),  can  be  calculated,  where,  from  Equa¬ 
tion  3.28,  the  initial  target  location  error,  y0,  is  equal  to  the  predicted  target  location 
error,  or 


Yo  t’npredO- 

The  resulting  combined  initial  state  vector  is 


(3.29) 
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Likewise,  the  resulting  combined  initial  covariance  matrix  is 


where  P^o  is  assigned  to  be 


(3.30) 


(3.31) 
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and  Pxyo,  Pyxo  and  Pyyo  are  calculated,  respectively,  as 


(3.32) 
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(3.33) 

(3.34) 

(3.35) 


The  influence  coefficient,  Gyx0,  defined  in  Equation  (2.75),  is  determined  from  the 
partial  derivatives  of  the  target  location  error  equation  solved  in  Section  3.4,  defined  as 


tnx  t'nxorb  T  tnxFOGM 
tny  tnyorb  T  tnyFOGM  • 


(3.36) 

(3.37) 
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Consider  the  homogeneous  nonlinear  differential  equation  defined  in  Equation  (2.16) 


x*(f)  =  f  [  x*(f), u*(f),  t  ].  (3.38) 

For  the  given  input  function,  uj)  (assumed  in  this  model  to  be  zero),  and  the  initial 
condition,  Xg,  the  nominal  solution  trajectory,  x*(f)  is  known  to  exist.  The  perturba¬ 
tions  in  the  initial  condition  are  denoted  as 

x*(t)  =  x*(t)  +  5x*(t),  (3.39) 

where  the  perturbation  state  vector  for  the  system  model  is  defined  as 
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From  Equations  (2.24)  and  (2.25),  the  transition  matrix  and  linear  perturbation  equa¬ 
tions,  respectively,  are  therefore 


$*(Mo)  =  eF*VAt  (3.41) 

5x*(t)  =  <3>*(f,f0)  5x*(t~)  +  w (t).  (3.42) 

From  Equations  (2.88)  through  (2.90)  and  Equations  (3.1)  through  (3.7),  the 
partial  derivative  dynamics  matrix  for  both  the  navigation  and  landmark  states.  F*(f), 
is  calculated  as 


F *(t)  = 


0100000  000 


|ik  o  0  0  0  000 

OT  x  OT  y  OT  z 


0001000  000 


^0^0^  0  0  000 

OT  x  OT  y  OT  z 


0000010  000 


l^o^o^oo  000 


dr 


oooooogooo 


oooooooffoo 


0000000  000 


0000000  000 


(3.43) 


F*(f)  can  be  further  solved  as 
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Note  that  the  dynamics  of  the  landmark  states,  tnx  and  tny,  are  zero,  since  they  are 
stochastic  projections  of  the  navigation  state  (see  Section  2.8).  From  Equations  (2.92) 
through  (2.94),  and  Equation  (3.13),  the  partial  derivative  pixel  projection  matrix, 
is  calculated  as 
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H*(i)  can  be  further  solved  as 
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where  H'  is  defined  as 
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The  resulting  estimates  of  the  navigation  and  landmark  state  will  be  presented 
in  the  next  chapter. 
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IV.  Results  and  Observations 


In  order  to  validate  the  mathematical  methods  presented  in  Chapter  III,  the  satellite- 
based  image-aided  navigation  system  algorithms  are  evaluated  using  simulation- 
based  analyses.  First,  the  development  of  the  simulation  will  be  described.  An  analy¬ 
sis  of  Monte  Carlo  results  for  a  satellite  system  in  both  low  Earth  and  Molniya  orbits 
will  follow. 

4-1  Simulation  Development 

The  performance  of  the  extended  Kalman  filter  built  in  Section  3.7  is  verified 
using  a  statistical  ensemble  of  sample  functions,  totaling  100  sample  functions  per 
Monte  Carlo  run.  Each  Monte  Carlo  run  has  a  20  minute  duration,  sampling  at 
1  Hz  (for  a  total  of  1200  data  points  per  sample).  For  each  sample,  an  entirely 
new  set  of  inertial  sensor  and  imaging  system  data  is  generated  using  the  respective 
system  dynamics  and  and  error  models  described  Chapter  III,  and  are  based  upon 
the  assigned  statistical  and  noise  parameters  listed  in  Table  3.6. 

For  each  20  minute  run,  the  EKF  is  implemented  for  three  distinct  profiles: 
the  first  without  any  image-aided  updates  to  the  target  location  error  estimate,  the 
second  using  only  the  image  matching  method  developed  in  Section  3.5  to  update 
the  target  location  error  estimate,  and  the  third  using  both  image  matching  and 
georeferenc-ing  methods  developed  in  Sections  3.5  and  3.6  to  update  the  target  location 
error  estimate.  For  convenience,  these  profiles  were  be  referred  to  as  “non-up dated”, 
“image-matched” ,  and  “image-corrected” ,  accordingly. 

In  each  profile,  initial  and  subsequent  target  location  error  predictions  are  made 
per  Equations  (3.28)  and  (3.29).  In  the  non- updated  profile,  these  predictions  are 
assigned  to  be  zero  mean  with  initial  uncertainty,  Pj).  In  the  image-matched  profile, 
subsequent  images  are  matched  to  the  first  captured  image  in  order  to  reduce  optical 
drift;  however,  without  knowledge  of  the  initial  target  location  error,  a  constant  target 
location  error  will  remain.  Finally,  in  the  image-corrected  profile,  the  initial  bias  of 
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the  image  matched  profile  is  corrected  through  georeferencing,  therefore  providing  the 
most  precise  target  geolocation  capability  of  the  three  profiles. 

In  order  to  demonstrate  system  performance  in  multiple  orbits,  two  Monte  Carlo 
scenarios  will  be  implemented.  The  first  scenario  contains  a  satellite  system  in  low 
Earth  orbit,  where  the  satellite  is  located  in  a  relatively  small  orbital  apogee  and  image 
resolution  is  high.  The  second  scenario  consist  of  a  system  in  a  Molniya  orbit,  where 
the  satellite  is  located  in  a  relatively  high  orbital  apogee  and  image  resolution  is  low. 
Orbital  and  imaging  system  parameters  are  listed  in  Tables  3.1  and  3.3,  respectively. 

For  simplicity,  the  pre-generated  high-resolution  and  low-resolntion  satellite  im¬ 
ages  of  the  greater  Cincinnati,  Ohio  area  are  used  to  represent  a  satellite  in  the  LEO 
and  MOL  orbits,  respectively.  This  template  is  depicted  in  Figure  3.4.  Whereas  these 
static  satellite  images  are  appropriate  for  testing  the  coupled  image/inertial  sensor  al¬ 
gorithm,  the  results  are  not  directly  comparable  to  the  performance  of  a  real  satellite 
imaging  system.  As  such,  imaging  issues  including  poor  sunlight  conditions,  target 
obstruction,  ego-motion  disparity  and  motion  blur  between  subsequent  images  are  not 
modeled  [25]. 

4-2  Low  Earth  Orbit  Simulation 

The  navigation  and  landmark  state  errors  of  the  imaging  satellite  are  simulated 
in  20  minute  durations  for  a  total  of  100  samples.  The  satellite  is  initialized  in  a  LEO 
orbit  at  an  elevation  of  598  km  (at  apogee)  above  the  Cincinnati,  Ohio  area,  with  an 
imaging  resolution  of  0.1  km/pixcl. 

The  position  and  velocity  errors  of  the  satellite,  in  the  e-frame  x-,  y-  and  z- 
directions,  are  shown  in  Figures  4.1  through  4.2,  respectfully.  As  expected  with  large 
acceleration-level  noise,  the  inertial  position  and  velocity  measurement  errors  accumu¬ 
late  over  time,  resulting  in  quickly  growing  position  and  velocity  error  uncertainties 
without  bound.  The  observed  drift  from  zero  mean  is  the  result  of  two  phenomenon. 
The  first  phenomenon  is  that  the  additive  error,  introduced  onto  the  initial  position 
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and  velocity  states  prior  to  propagation,  results  in  a  slight  deviation  of  the  estimated 
orbit  from  true  orbit  over  time,  producing  a  non-zero  position  and  velocity  error 
mean.  The  second  phenomenon  is  due  to  an  insufficient  number  of  Monte  Carlo  sam¬ 
ples,  resulting  in  a  perceived  “bias”  of  each  position  and  velocity  error  state  away 
from  the  true  zero  mean.  Although  not  confirmed  by  this  research,  it  is  hypothesized 
that  by  running  tens  of  thousands  of  samples,  the  position  and  velocity  errors  would 
begin  to  approach  their  expected  zero  means.  Note  the  variations  of  position  and  ve¬ 
locity  errors  between  non-updated,  image-matched  and  image-corrected  profiles  were 
minimal;  therefore,  only  the  non-updated  profile  is  depicted. 

The  attitude  errors  of  the  imaging  system,  time-correlated  first-order  Gauss- 
Markov  processes  described  in  Sections  3.2.1  and  3.3,  are  shown  in  Figure  4.3.  Being 
first-order  Gauss-Markov  processes,  the  state  estimates  are  zero-mean  with  growing 
uncertainties,  reaching  a  steady  state  of  approximately  0.17  degs  at  time  constant,  r, 
of  10  minutes.  Note  the  variations  of  attitude  errors  between  non-updated,  image- 
matched  and  image-corrected  profiles  were  minimal;  therefore,  only  the  non-updated 
profile  is  depicted. 

The  n-frarne  target  location  errors  are  estimated  with  respect  to  the  non- 
updated,  image-matched  and  image-corrected  profiles  as  shown  in  Figures  4.4  through 
4.6,  respectfully.  As  a  stochastic  projection  of  the  navigation  states,  the  behavior  of 
the  non-updated  target  location  errors  in  Figure  4.4  correspond  to  the  growing  posi¬ 
tion  and  velocity  errors  and  of  the  first-order  Gauss-Markov  attitude  errors.  As  such, 
a  zero-mean  target  location  error  is  observed  with  a  large,  unbounded  uncertainty 
growing  over  the  20  min  duration.  Recalling  Equation  (3.37),  the  target  location  er¬ 
ror  due  to  satellite  orbital  drift  is  dominate  in  LEO  orbit,  primarily  the  result  to  faster 
satellite  trajectories,  larger  off-nadir  angles  and  larger  atmospheric  effect  considera¬ 
tions,  causing  a  larger  initial  uncertainty  in  a  particular  direction  (in  this  case,  the 
x-direction).  The  non-updated  profile  can  be  visualized  with  respect  to  the  satellite 
image  in  Figure  4.7. 
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The  image-matched  target  location  errors  in  Figure  4.5  are  also  zero-mean;  how¬ 
ever.  it  can  be  seen  that  the  time-correlated  drift  with  respect  to  the  attitude  states  is 
suppressed  by  the  image  matching  algorithm,  resulting  in  steady  uncertainties.  This 
uncertainty  remains  high,  however,  since  navigation  state  uncertainties  are  large  and 
knowledge  of  the  initial  target  location  error  is  not  available.  Again,  the  satellite  or¬ 
bital  drift  clue  to  large  statistical  and  noise  parameters  dominates  the  target  location 
error,  potentially  causing  a  larger  uncertainty  in  a  particular  direction  (in  this  case, 
the  x-direction).  The  image-matched  profile  can  be  visualized  with  respect  to  the 
satellite  image  in  Figure  4.8. 

In  Figure  4.6,  the  image-corrected  target  location  error  remains  small  in  both 
x-  and  y-directions.  This  error  is  seen  to  improve  by  roughly  an  order  of  magnitude 
compared  to  that  of  the  previous  two  profiles.  The  image-corrected  profile  takes 
advantage  of  both  the  image-matching  algorithm  (which  corrects  the  time-correlated 
drift)  as  well  as  the  georeferencing  algorithm  (which  corrects  the  initial  target  location 
error).  The  result  is  a  greatly  reduced  target  location  error  throughout  the  entire  20 
minute  duration,  corrupted  only  by  minor  measurement  noise.  The  image-corrected 
profile  can  be  visualized  with  respect  to  the  satellite  image  in  Figure  4.9. 

Finally,  the  root-sum-squared  (RSS)  errors  of  the  target  location  error  are  ana¬ 
lyzed  in  order  to  provide  a  more  direct  comparison  of  the  simulated  satellite  system’s 
performance  with  respect  to  the  three  profiles.  The  RRS  errors  comparing  the  tar¬ 
get  location  errors  in  these  three  cases  are  shown  in  Figure  4.10.  Over  the  entire  20 
minute  simulation,  it  can  be  clearly  seen  that  the  image-corrected  profile  improves 
the  system  performance  by  an  order  of  magnitude  over  that  of  the  non-updated  and 
image-matched  profiles. 
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Figure  4.1:  Simulated  100-run  Monte  Carlo  satellite  position  error  results  in  LEO, 
without  image  updates.  The  position  error  sample  functions  are  indicated 
by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard  deviation 
are  indicated  by  the  green  and  red  solid  lines,  respectively.  The  large 
acceleration-level  noise  results  in  quickly  growing  position  uncertainties. 
The  observed  drift  is  most  likely  due  to  both  the  initial  introduction  of 
acceleration  noise  as  well  as  the  low  number  of  Monte  Carlo  samples. 
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Figure  4.2:  Simulated  100-run  Monte  Carlo  satellite  velocity  error  results  in  LEO, 
without  image  updates.  The  velocity  error  sample  functions  are  indicated 
by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard  deviation 
are  indicated  by  the  green  and  red  solid  lines,  respectively.  The  large 
acceleration-level  noise  results  in  quickly  growing  velocity  uncertainties. 
The  observed  drift  is  most  likely  due  to  both  the  initial  introduction  of 
acceleration  noise  as  well  as  the  low  number  of  Monte  Carlo  samples. 
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Figure  4.3:  Simulated  100-run  Monte  Carlo  satellite  attitude  error  results  in  LEO, 
without  image  updates.  The  attitude  error  sample  functions  are  indicated 
by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard  deviation 
are  indicated  by  the  green  and  red  solid  lines,  respectively.  The  attitude 
errors  are  zero-mean  with  growing  uncertainties,  reaching  a  steady  state 
of  approximately  0.17  degs  at  a  time  constant,  r,  of  10  minutes. 
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Figure  4.4:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  LEO,  with¬ 
out  image  updates.  The  target  location  error  sample  functions  are  indi¬ 
cated  by  blue  dotted  lines.  The  ensemble  mean  and  standard  deviation  are 
indicated  by  the  green  and  red  solid  lines,  respectively.  A  zero-mean  tar¬ 
get  location  error  is  observed  with  a  large,  unbounded  uncertainty  growing 
over  the  20  min  duration. 
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Figure  4.5:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  LEO,  with 
image  matching.  The  target  location  error  sample  functions  are  indicated 
by  blue  dotted  lines.  The  ensemble  mean  and  standard  deviation  are  in¬ 
dicated  by  the  green  and  red  solid  lines,  respectively.  The  time-correlated 
drift  with  respect  to  the  attitude  errors  is  suppressed  by  the  image  match¬ 
ing  algorithm,  resulting  in  large  steady  uncertainties. 
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Figure  4.6:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  LEO,  with 
full  image  correction.  The  target  location  error  sample  functions  are  in¬ 
dicated  by  blue  dotted  lines.  The  ensemble  mean  and  standard  deviation 
are  indicated  by  the  green  and  red  solid  lines,  respectively.  Both  time- 
correlated  drift  and  initial  target  location  error  are  corrected  in  this  profile, 
resulting  in  a  reduced  target  location  error  by  an  order  of  magnitude. 
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Figure  4.7:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  minutes 
in  LEO,  without  image  updates.  Over  time,  the  target  location  error  is 
seen  to  grow  without  bound. 
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Figure  4.8:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  minutes 
in  LEO,  with  image  matching.  Although  the  initial  target  location  error 
is  large,  the  error  does  not  significantly  drift  over  time. 
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Figure  4.9:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  minutes 
in  LEO,  with  full  image  correction.  Both  the  image- matching  algorithm 
(which  corrects  the  time-correlated  drift)  and  georeferencing  algorithm 
(which  corrects  the  initial  target  location  error)  are  utilized. 
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Figure  4.10:  Simulated  100-run  Monte  Carlo  root-sum-squared  (RSS)  target  location 
error  results  in  LEO,  comparing  the  three  image-aided  profiles.  It  can 
be  clearly  seen  that  the  image-corrected  profile  improves  the  system  per¬ 
formance  by  an  order  of  magnitude  over  that  of  the  non-updated  and 
image-matched  profiles. 
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4-3  High  Earth  Orbit  Simulation 

For  a  satellite  system  in  a  MOL  orbit,  an  identical  simulation  to  that  of  the 
LEO  scenario  is  conducted  to  determine  whether  the  satellite  system’s  elevation  and 
image  resolution  would  vastly  affect  the  results.  For  comparison,  the  navigation  and 
landmark  state  errors  of  the  imaging  satellite  are  again  simulated  in  20  min  durations 
for  a  total  of  100  samples.  The  satellite  is  initialized  in  a  MOL  orbit  at  an  elevation 
of  38,900  km  (at  apogee)  above  the  Cincinnati,  Ohio  area,  with  an  imaging  resolution 
of  1  km/pixel. 

The  e-frame  x-,y-  and  z-directional  position  and  velocity  errors  of  the  satellite 
are  shown  in  Figures  4.11  through  4.12,  respectfully.  With  the  smaller  acceleration- 
level  noise  at  the  higher  orbit,  the  inertial  position  and  velocity  measurement  errors 
accumulate  more  slowly  over  time,  resulting  in  slowly  growing  position  and  velocity 
error  uncertainties  without  bound.  As  with  the  LEO  orbit  position  and  navigation 
error  states,  the  observed  drift  is  due  to  both  the  initial  introduction  of  acceleration 
noise  as  well  as  the  low  number  of  Monte  Carlo  samples.  Again,  the  variations  of  po¬ 
sition  and  velocity  errors  between  non-updated,  image-matched  and  image-corrected 
profiles  were  minimal;  therefore,  only  the  non-updated  profile  is  depicted. 

The  time-correlated  first-order  Gauss-Markov  process  attitude  errors  of  the 
imaging  system  are  shown  in  Figure  4.13.  Being  first-order  Gauss-Markov  processes, 
the  state  estimates  are  zero-mean  with  growing  uncertainties,  reaching  a  steady  state 
of  approximately  0.017  degs  at  time  constant,  r,  of  10  minutes.  Again,  the  variations 
of  attitude  errors  between  non-updated,  image-matched  and  image-corrected  profiles 
were  minimal;  therefore,  only  the  non-updated  profile  is  depicted. 

The  n-frame  target  location  errors  are  estimated  with  respect  to  the  non- 
updated,  image- matched  and  image-corrected  profiles  as  shown  in  Figures  4.14  through 
4.16,  respectfully.  As  with  the  lower  orbit  scenario,  the  behavior  of  the  non-updated 
target  location  errors  in  Figure  4.14  correspond  to  the  growing  position,  velocity,  and 
attitude  errors.  As  such,  a  zero-mean  target  location  error  is  observed  with  a  large, 
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unbounded  uncertainty  growing  over  the  20  min  duration.  The  target  location  error  in 
the  MOL  scenario  is  not  dominated  by  either  satellite  orbital  drift  or  time-correlated 
optical  drift,  resulting  in  a  large  and  unsteady  initial  uncertainty  in  both  x-  or  y- 
directions.  The  non-updated  profile  can  be  visualized  with  respect  to  the  satellite 
image  in  Figure  4.17. 

The  image-matched  target  location  errors  in  Figure  4.15  are  also  zero-mean; 
however,  like  in  the  LEO  scenario,  the  time-correlated  attitude  drift  is  suppressed  by 
the  image  matching  algorithm.  This  results  in  steady  and  high  uncertainties  dne  to 
large  navigation  state  uncertainties  and  unknown  initial  target  location  error.  In  the 
MOL  orbit,  neither  orbital  drift  or  optical  drift  dominate  the  target  location  error. 
The  image-matched  profile  can  be  visualized  with  respect  to  the  satellite  image  in 
Figure  4.18. 

In  Figure  4.16,  the  image-corrected  target  location  error  remains  small  in  both 
x-  and  y-directions.  Like  the  LEO  scenario,  this  error  is  seen  to  improve  by  roughly 
an  order  of  magnitude  compared  to  that  of  the  previous  two  profiles,  taking  advantage 
of  both  the  image-matching  and  georeferencing  algorithms.  The  result  is  a  greatly 
reduced  target  location  error  throughout  the  entire  20  minute  duration,  corrupted 
only  by  minor  measurement  noise.  The  image-corrected  profile  can  be  visualized  with 
respect  to  the  satellite  image  in  Figure  4.19. 

Finally,  the  RSS  errors  of  the  target  location  error  are  analyzed  with  respect  to 
the  three  profiles.  The  RRS  errors  comparing  the  target  location  errors  in  these  three 
cases  are  shown  in  Figure  4.20.  Very  similar  to  the  LEO  scenario,  it  can  be  clearly 
seen  that  the  image-corrected  profile  improves  the  system  performance  by  an  order 
of  magnitude  over  the  entire  20  minute  duration. 

4-4  Simulation  Comparisons 

In  comparison  of  the  LEO  and  MOL  scenarios,  the  image-corrected  profile  in 
both  cases  has  been  shown  to  provide  highly  accurate  target  tracking  results  over  a 
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20  minute  duration.  The  larger  statistical  and  noise  parameters  of  the  LEO  orbit, 
resulting  in  respectively  large  navigation  uncertainties,  still  do  not  drastically  affect 
target  tracking  performance  in  the  image-matched  or  image-corrected  profiles.  Like¬ 
wise,  in  both  the  high-resolution  LEO  system  and  low-resolution  MOL  system,  it  is 
shown  that  subsequent  satellite  images  can  be  effectively  matched  and  minimal  initial 
target  location  error  can  be  provided  through  georeferencing,  in  order  to  accurately 
predict  and  track  the  target  location  in  either  scenario. 
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Figure  4.11:  Simulated  100-run  Monte  Carlo  satellite  position  error  results  in  MOL, 
without  image  updates.  The  position  error  sample  functions  are  indi¬ 
cated  by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard 
deviation  are  indicated  by  the  green  and  red  solid  lines,  respectively.  The 
smaller  acceleration- level  noise  results  in  slowly  growing  position  uncer¬ 
tainties.  The  observed  minor  drift  is  most  likely  due  to  both  the  initial 
introduction  of  acceleration  noise  as  well  as  the  low  number  of  Monte 
Carlo  samples. 
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Figure  4.12:  Simulated  100-run  Monte  Carlo  satellite  velocity  error  results  in  MOL, 
without  image  updates.  The  velocity  error  sample  functions  are  indicated 
by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard  deviation 
are  indicated  by  the  green  and  red  solid  lines,  respectively.  The  smaller 
acceleration-level  noise  results  in  very  minor  growth  of  velocity  error 
uncertainties. 
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Figure  4.13:  Simulated  100-run  Monte  Carlo  satellite  attitude  error  results  in  MOL, 
without  image  updates.  The  attitude  error  sample  functions  are  indi¬ 
cated  by  blue  dotted  lines.  The  ensemble  mean  and  3-sigma  standard 
deviation  are  indicated  by  the  green  and  red  solid  lines,  respectively. 
The  attitude  errors  are  zero-mean  with  growing  uncertainties,  reaching 
a  steady  state  of  approximately  0.017  degs  at  a  time  constant,  r,  of  10 
minutes. 
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Figure  4.14:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  MOL, 
without  image  updates.  The  target  location  error  sample  functions  are 
indicated  by  blue  dotted  lines.  The  ensemble  mean  and  standard  de¬ 
viation  are  indicated  by  the  green  and  red  solid  lines,  respectively.  A 
zero-mean  target  location  error  is  observed  with  a  large,  unsteady  un¬ 
certainties  over  the  20  min  duration. 
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Figure  4.15:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  MOL, 
with  image  matching.  The  target  location  error  sample  functions  are 
indicated  by  blue  dotted  lines.  The  ensemble  mean  and  standard  devi¬ 
ation  are  indicated  by  the  green  and  red  solid  lines,  respectively.  The 
time-correlated  drift  with  respect  to  the  attitude  errors  is  suppressed  by 
the  image  matching  algorithm,  resulting  in  large  steady  uncertainties. 
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Figure  4.16:  Simulated  100-run  Monte  Carlo  target  location  error  results  in  MOL, 
with  full  image  correction.  The  target  location  error  sample  functions 
are  indicated  by  blue  dotted  lines.  The  ensemble  mean  and  standard 
deviation  are  indicated  by  the  green  and  red  solid  lines,  respectively. 
Both  time-correlated  drift  and  initial  target  location  error  are  corrected 
in  this  profile,  resulting  in  a  reduced  target  location  error  by  an  order  of 
magnitude. 
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Figure  4.17:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  min¬ 
utes  in  MOL,  without  image  updates.  Over  time,  the  target  location 
error  is  seen  to  very  slowly  grow  without  bound. 


Figure  4.18:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  min¬ 
utes  in  MOL,  with  image  matching.  Although  the  initial  target  location 
error  is  large,  the  error  does  not  significantly  drift  over  time. 
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Figure  4.19:  Mapped  100-run  Monte  Carlo  target  location  error  results  after  20  min¬ 
utes  in  MOL,  with  full  image  correction.  Both  the  image-matching  al¬ 
gorithm  (which  corrects  the  time-correlated  drift)  and  georeferencing 
algorithm  (which  corrects  the  initial  target  location  error)  are  utilized. 


100 


RSS  T arget  Loc  Err  (pix) 


Figure  4.20:  Simulated  100-run  Monte  Carlo  root-sum-squared  (RSS)  target  location 
error  results  in  MOL,  comparing  the  three  image-aided  profiles.  The 
image-corrected  profile  improves  the  system  performance  by  an  order  of 
magnitude  over  that  of  the  non-updated  and  image-matched  profiles. 
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V.  Conclusions 


This  thesis  introduces  the  concept  of  fusing  the  imaging  and  inertial  sensors  of 
a  satellite  observation  system  to  accurately  and  autonomously  geolocate  com¬ 
mercial  or  military  ground  targets  of  interest.  In  this  chapter,  conclusions  regarding 
the  image-aided  satellite  system  simulation  are  presented,  and  potential  focus  areas 
of  future  research  are  addressed. 

5. 1  Summary 

As  mentioned  in  Chapter  I,  the  goal  of  this  research  was  to  develop  a  low-cost, 
low- weight,  highly- accurate  image-aided  inertial  satellite  navigation  system  without 
the  need  of  human  interaction  or  dependency  on  external  navigation  system  sources. 
This  section  summarizes  the  implementation,  results  and  observations  of  this  system 
design. 

In  Chapter  III,  the  orbital  modeling  parameters  listed  in  Section  3.1  were  used  to 
define  both  a  nominal  low  Earth  and  Molniya  orbit  of  an  imaging  satellite.  The  satel¬ 
lite  system  parameters  and  dynamics  defined  in  Section  3.2  described  the  trajectory  of 
the  vehicle  and  the  functionality  of  the  onboard  imaging  system.  Next,  the  navigation 
and  landmark  noise  error  parameters  assigned  in  Section  3.3  were  introduced  into  the 
nominal  state  initial  conditions  and  system  dynamics  to  produce  simulated  truth  data; 
specifically,  the  generated  ground  error  as  described  in  Section  3.4.  Image  matching 
and  georeferencing  techniques  presented  in  Sections  3.5  and  3.6,  respectively,  aided 
in  the  tracking  of  the  ground  target  of  interest  by  detecting  unique  features  of  the 
target  between  subsequent  images  and  correcting  for  any  initial  target  location  error. 
Finally,  the  extended  Kalman  filter  described  in  Section  3.7  was  implemented  in  order 
to  estimate  of  the  navigation  and  landmark  states  and  determine  the  errors  between 
the  estimated  and  truth  data  over  time. 

In  Chapter  IV,  the  satellite-based  image-aided  navigation  system  algorithms  de¬ 
fined  in  Chapter  III  were  evaluated  using  Monte  Carlo  simulation-based  analyses.  The 
performance  of  the  extended  Kalman  filter  was  verified  using  a  20  minute  statistical 
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ensemble  of  100  independent  sample  functions,  sampled  at  1  Hz,  and  were  analyzed  in 
three  image-aided  profiles:  non-updated,  image-matched  and  image-corrected.  Two 
separate  Monte  Carlo  scenarios  are  implemented:  the  first  was  a  satellite  system  in  a 
low  Earth  orbit  (with  high  image  resolution)  and  the  second  system  in  a  high  Earth 
orbit  (with  low  image  resolution).  Respective  high  and  low  resolution  images  of  the 
Cincinnati,  Ohio  area  were  used  to  represent  the  satellite’s  field  of  view  and  identify  a 
unique  target  of  interest.  In  both  LEO  and  MOL  orbit  scenarios  (Sections  4.2  and  4.3, 
respectively),  it  was  observed  that  the  full  image-aided  target  location  prediction  of 
the  image-corrected  profile  minimized  the  target  location  error  by  an  order  of  mag¬ 
nitude  throughout  the  20  minute  simulation.  Additionally,  it  was  shown  that  the 
satellite  system’s  trajectory  and  image  resolution  capability  did  not  drastically  affect 
the  performance  of  the  image-aided  satellite  system  in  either  scenario,  verifying  that 
this  model  is  suitable  for  both  high  or  low  Earth  orbit  assignments  with  corresponding 
image  resolution  requirements. 

5.2  Conclusions 

Based  upon  the  methodology,  results  and  observations  summarized  in  Sec¬ 
tion  5.1,  final  conclusions  of  the  image-aided  satellite  inertial  navigation  system  can 
be  made. 

A  significant  motivator  defined  in  Chapter  I  expressed  the  need  for  a  low-cost, 
low-weight,  highly-accurate  satellite  imaging  system.  Growing  expenses  associated 
with  overall  satellite  design  and  space  launch  capabilities  demand  that  space  vehicles 
be  as  efficient  and  light  weight  as  possible.  Since  the  developed  image-aided  algo¬ 
rithms  defined  in  Chapter  III  utilize  only  pre-existing  image  and  inertial  sensors,  and 
no  additional  vehicle  or  ground  tracking  hardware  is  implemented,  it  can  be  reason¬ 
ably  concluded  that  vehicle  cost  and  weight  will  not  be  drastically  affected  by  this 
design.  The  accuracy  of  the  fully  corrected  image-aided  model,  as  discussed  in  detail 
in  Chapter  IV,  provides  appreciable  evidence  that  minimal  target  location  error  is 
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available  for  the  first  20  minutes  of  target  tracking.  It  is  assumed  that  this  period  of 
time  is  sufficient  for  most  target  tracking  applications. 

Since  this  integrated  image/inertial  Earth  observation  system  actually  estimates 
the  extended  Kalman  filter  navigation  states  of  the  satellite,  these  estimates  can  be 
fed  back  into  the  satellite  in  order  to  improve  the  combined  satellite  navigation/target 
tracking  solution.  This  is  particularly  beneficial  to  navigation  systems  requiring  rou¬ 
tine  navigation  calibrations  since  these  calibrations  are  accomplished  simultaneously 
with  the  target  tracking  functionality. 

Another  motivator  for  this  research  was  to  build  a  system  that  is  fully  au¬ 
tonomous  and  independent  of  external  navigation  reference  sources,  such  as  GPS  or 
continuous  post-capture  georeferencing.  The  image  matching  SSD  algorithm  is  shown 
to  be  fully  autonomous,  and  therefore  satisfies  this  requirement.  Likewise,  GPS  or 
similar  external  reference  sources  are  never  implemented  in  this  model  and  is  therefore 
independent  of  such  systems.  The  georeferencing  algorithm  used  to  initially  correct 
the  target  location  prediction  could  be  viewed  as  a  possible  infringement;  however,  in 
this  application,  it  is  assumed  that  a  pre-generated  georeferenced  map  of  the  area  of 
interest  is  available,  with  a  sufficient  number  of  control  points  to  accurately  correct 
any  initial  target  location  errors.  Additionally,  since  an  image  matching  algorithm 
is  used  in  conjunction  with  this  georeferencing  technique,  the  target  need  only  be 
georeferenced  once  in  order  to  provide  significant  target  tracking  improvement.  As 
a  comparison,  if  the  target  were  to  instead  be  georeferenced  once  every  sample  (in 
the  case  of  no  available  image- matching  capabilities),  sluggish  georeferencing  response 
time  would  potentially  outweigh  any  target  location  error  correction.  It  was  shown 
that  one-time  geolocation  early  in  the  simulation  provided  considerable  target  track¬ 
ing  improvement  throughout  the  entire  simulation;  therefore,  the  response  time  of  one 
georeferencing  sampling  is  not  of  major  concern. 
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Although  further  advances  could  be  made  to  further  optimize  the  performance 
of  the  target  tracking  system,  these  results  show  promise  that  the  development  of  a 
robust  image-aided  satellite  system  is  worthwhile. 

5.3  Future  Work 

This  research  represents  a  preliminary  analysis  in  the  fusion  of  image/inertial 
sensors  of  Earth-observation  satellites  for  precise  geolocation  of  ground  targets.  A 
number  of  recommendations  for  future  research  can  be  made,  further  exploring  this 
space-based  image-aided  navigation  system.  They  are  listed  as  follows: 

•  Simulation  Using  Real  Data:  While  the  pre-generated  satellite  images  are  ap¬ 
propriate  for  testing  the  coupled  image/inertial  sensor  algorithm,  the  results  are 
not  directly  comparable  to  the  performance  of  a  real  satellite  imaging  system. 
A  logical  next  step  would  be  to  test  a  series  of  authentic,  subsequent  satellite 
images  of  a  feature-rich  ground  target  in  order  to  provide  further  verification  of 
system  performance. 

•  Use  of  a  Standard  Earth  Coordinate  Frame:  This  research  uses  a  simplistic, 
spherical  Earth  model  as  a  reference  frame  and  assumes  the  target  of  interest  is 
located  at  a  fixed  elevation  (i.e.,  at  sea  level).  A  more  sophisticated  reference 
coordinate  system,  such  as  the  World  Geodetic  System  1984  (WGS  84),  would 
provide  more  precise  ground  tracking  capabilities  for  real-world  image-aided 
satellite  navigation  systems. 

•  Introduction  of  Imaging  Anomalies:  Imaging  issues  such  as  poor  sunlight  con¬ 
ditions,  partial  target  obstruction,  binocular  disparity,  motion  blur,  and  affine 
transformations  are  not  identified  in  this  research.  It  is  recommended  that  a 
number  of  these  considerations  be  modeled  to  provide  further  system  credibility 
in  the  presence  of  real-world  imaging  anomalies. 
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